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Abstract 

As  navigation  systems  continue  to  improve  in  performance  and  features,  the  Air  Eorce 
must  develop  better  Navigation  Reference  Systems  (NRS)  to  keep  pace  with  technology. 
Specifically,  with  the  advent  of  enhanced,  integrated  Global  Positioning  System  (GPS) 
and  Inertial  Navigation  System  (INS)  navigators,  emphasis  is  placed  on  the  measuring 
performance  in  the  presence  of  GPS  jamming.  To  meet  these  needs,  a  new  NRS,  dubbed 
the  Sub-Meter  Accuracy  Reference  System  (SARS),  is  being  developed  by  the  746th  Test 
Squadron,  Holloman  AEB,  New  Mexico.  SARS  uses  a  unique,  inverted  GPS  pseudolite 
positioning  system  to  determine  a  reference  trajectory.  This  research  investigates  two  post¬ 
processing  methods  of  determining  velocity  from  a  discrete  position  data  at  a  constant 
data  rate.  The  first  method  employs  numerical  differentiation  along  with  digital  filters  to 
provide  noise  reduction.  The  second  method  uses  kinematic  model-based  Kalman  filtering 
and  smoothing  to  determine  the  reference  velocity. 
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VELOCITY  DETERMINATION  FOR  AN  INVERTED  PSEUDOLITE 
NAVIGATION  REFERENCE  SYSTEM 

/.  Introduction 

I. 1  Background 

The  746th  Guidance  Test  Squadron,  Holloman  AFB,  NM,  is  the  principal  test  or¬ 
ganization  for  developmental  testing  of  military  aircraft  navigation  systems  [39].  As  part 
of  the  746th,  the  Central  Inertial  Guidance  Test  Facility  (CIGTF)  maintains  a  Navigation 
Reference  System  (NRS)  used  in  the  test  and  evaluation  of  navigation  systems. 

From  1965  to  the  present,  many  different  technologies  have  been  used  in  the  de¬ 
velopment  of  an  NRS.  These  include  radar  tracking,  high-precision  ground-based  camera 
tracking  and  the  development  of  an  aircraft  transponder /ground  receiver  system.  In  1975, 
with  the  advent  of  microprocessor  technology,  advances  in  mathematics  (Kalman  filtering) 
and  Inertial  Navigation  Systems  (INS),  the  “modern-era”  NRS  was  conceived.  The  Com¬ 
pletely  Integrated  Reference  Instrumentation  System  (CIRIS)  combined  an  INS,  baromet¬ 
ric  altimeter  and  a  Range/ Range- Rate  System  (RRS)  of  ground  transponder/interrogators. 
Over  the  years,  CIRIS  was  updated  with  new  inertial  navigators  as  well  as  Global  Posi¬ 
tioning  System  (GPS)  receivers.  These  upgrades  continually  improved  the  overall  accuracy 
of  the  system. 

In  1993,  development  began  on  a  replacement  for  CIRIS  called  the  CIGTF  High 
Accuracy  Post-Processing  Reference  System  (CHAPS).  CHAPS  combined  an  INS,  Differ¬ 
ential  GPS  and  the  RRS  transponder  measurements  in  a  Kalman  filter,  in  a  manner  similar 
to  CIRIS.  CHAPS,  however,  had  to  be  significantly  more  accurate  than  CIRIS  in  order  to 
effectively  test  GPS  systems,  inertial  systems  and  GPS-aided  inertial  navigation  systems. 
As  an  NRS,  CHAPS  has  shown  performance  in  position  accuracy  of  1.2  m  and  velocity 
accuracy  of  0.03  m/s  [39].  In  early  1995,  an  upgraded  version  of  CHAPS,  called  CHAPS 

II,  began  validation  testing  at  CIGTF.  CHAPS  II  is  a  more  compact  version  of  CHAPS 
with  an  enhanced  GPS  receiver  capable  of  carrier-phase  positioning. 
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To  be  an  effective  reference,  the  NRS  must  demonstrate  an  accuracy  at  least  one  or¬ 
der  of  magnitude  better  than  the  system  under  test  [2].  Advances  in  navigation  technology 
have  made  maintaining  the  accuracy  advantage  of  the  NRS  difficult.  Today,  technology 
continues  to  improve  the  accuracy  of  navigation  systems.  Embedded  GPS/INS  systems, 
designed  for  enhanced  performance  in  routine  flight,  high  dynamics,  antenna  shaded  and 
jamming  environments,  are  coming  to  market.  One  such  system,  the  Honeywell  H-764G, 
has  demonstrated  position  accuracy  of  0.8  nmi/hr  circular  error  probable  (CEP)  and  ve¬ 
locity  accuracy  of  0.06  mjs  root  mean  square  (RMS),  using  P-code  GPS  [12].  While  a 
pre-production  version  of  the  H-764G  has  been  evaluated  at  CIGTF  [30],  the  increased  po¬ 
sition  and  velocity  accuracy  of  such  systems,  along  with  claimed  anti-jamming  capabilities 
necessitates  a  performance  enhancement  of  the  current  NRS,  CHAPS.  A  comparison  of 
recent  reference  systems  developed  at  CIGTF  is  given  in  Table  1.1.  A  summary  of  earlier 
reference  system  development  can  be  found  in  [41]. 


Table  1.1  Comparison  of  Modern  Navigation  Reference  Systems 


Reference 

System 

Measurements 

Position 

Accuracy 

Velocity 

Accuracy 

CIRIS  [41] 

INS 

RRS 

GPS 

4.0  m 

0.03  m/s 

CHAPS  [36] 

INS 

RRS 

Differential  GPS 

1.2  m 

0.03  m/s 

CHAPS  II  [36] 

INS 

RRS 

Carrier-Phase  GPS 

0.1  m 

0.03  m/s 

SARS 

Pseudolites 

0.1  m 

Holloman  High  Speed 
Test  Track  [36] 

Instrumented  track 

0.0008  m 

0.0004  m/s 

To  meet  the  requirements  for  the  NRS,  CIGTF  is  developing  a  Sub-meter  Accuracy 
Reference  System  (SARS).  SARS  is  being  designed  for  accuracies  on  the  order  of  0.1  m 
in  position,  and  0.005  m/s  in  velocity  [13].  SARS  is  expected  to  operate  in  the  most 
dynamic  environments  experienced  by  manned  fighter  aircraft,  up  to  9  gees  and  10  m/s^ 


^Desired  accuracies. 


jerk.  CHAPS,  on  the  other-hand,  was  designed  to  fly  on  cargo  or  large  bomber  aircraft  [39]. 
CHAPS  II,  while  one-half  the  size  of  CHAPS,  is  unable  to  fly  aboard  operational  fighter 
aircraft  due  to  it’s  size.  Additionally,  SAKS  intends  to  operate  without  degradation  under 
narrow-band  GPS  jamming  conditions. 

To  accomplish  these  objectives,  SARS  departs  radically  from  the  NRS  designs  of 
CIRIS  and  CHAPS^.  SARS  will  depend  solely  on  GPS-style  measurements  derived  from 
the  carrier-phase  observable  [34].  Instead  of  using  GPS  receivers  on-board  the  aircraft,  the 
aircraft  will  be  fitted  with  a  GPS  pseudolite  transmitter,  which  is  explained  in  Section  1.3.3. 
This  dramatically  decreases  the  power  and  space  requirements  of  the  reference  system 
aboard  the  test  aircraft,  allowing  for  a  broader  range  of  aircraft  to  be  tested.  Furthermore, 
the  use  of  customizable  GPS  pseudolites  allows  the  user  to  change  the  carrier  frequency 
of  the  transmitters.  This  may  allow  testing  under  narrow-band  GPS  jamming  conditions, 
without  adversely  effecting  the  reference  system. 

1.2  Problem 

In  previous  NRS  designs,  an  error-state  Kalman  filter  combines  the  available  mea¬ 
surements  from  navigation  aids  such  as  GPS,  to  estimate  INS  errors.  The  INS  measures 
specific  force.  Acceleration  of  the  aircraft  is  obtained  by  subtracting  the  gravity  com¬ 
ponent  from  measured  specific  force.  Integrating  once  yields  velocity;  integrating  twice 
yields  position.  SARS,  on  the  other-hand,  will  obtain  position  information  from  GPS  data 
only,  via  the  GPS  observation  equations.  One  might  approximate  velocity  by  numerical 
differentiation,  however  at  the  expense  that  noise  present  in  the  position  measurements  is 
amplified  by  the  differentiation  procedure. 

In  this  effort  we  will  compare  two  different  approaches  to  determining  the  velocity 
of  an  aircraft  along  a  trajectory  given  GPS  measurements  of  position.  The  first  approach 
uses  a  combination  of  digital  filtering  techniques  and  numerical  methods.  The  second 
approach  uses  simple  kinematic  models  in  a  Kalman  filter  to  provide  estimates  of  velocity. 
We  will  compare  the  performance  of  these  two  methods  against  the  desired  performance 

^Conceptually,  the  SARS  design  is  similar  to  the  transponder  portion  of  the  RRS,  which  is  part  of 
CHAPS  and  CIRIS.  The  RRS  was  originally  developed  for  CIGTF  by  the  Cubic  Corporation,  in  1973  [41]. 
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specification  of  0.005  m/s  and  determine  the  sensitivity  of  each  algorithm  to  different 
parameters  such  as  measurement  rates  and  noise  characteristics.  A  simulated  environment 
will  be  used  for  this  research. 

1.3  Summary  of  Current  Knowledge 

1.3.1  Terminology.  Due  to  redundant  and  sometimes  inconsistent  terms  found 
in  the  open  literature,  it  is  necessary  to  untangle  the  terminology  used  in  discussing  GPS 
measurements.^  Three  “raw”  measurements  are  obtainable  in  high  performance  GPS 
receivers:  pseudorange,  deltarange  and  the  carrier-phase  observable.  Pseudorange,  also 
known  as  the  code  phase  observable,  is  the  traditional  GPS  measurement  used  to  deter¬ 
mine  position.  Using  a  code  correlation  technique,  pseudorange  is  a  measure  of  the  time 
of  transmission  between  the  GPS  satellite  and  receiver,  expressed  in  terms  of  distance.  A 
minimum  of  four  unique  pseudorange  measurements  are  required  to  compute  a  3-D  solution 
for  user  position  and  time. 

The  carrier-phase  observable  is  also  employed  to  determine  position.  Carrier-phase 
is  a  measurement  of  the  number  of  integer  and  fractional  cycles  of  the  GPS  carrier  signal 
between  the  transmitter  and  receiver,  expressed  in  terms  of  distance.  Using  carrier-phase 
to  resolve  position  with  centimeter  accuracy  was  first  developed  for  static  applications, 
such  as  surveying.  Recently,  the  ability  to  resolve  the  carrier-phase  ambiguity  problem  in 
real-time  has  been  demonstrated  [14,18,19],  leading  the  way  for  the  use  of  carrier-phase 
in  dynamic  environments. 

Deltarange,  also  known  as  Doppler,  differs  from  pseudorange  and  carrier-phase,  as  it 
is  a  measure  of  range-rate  rather  than  position.  Deltarange  approximates  the  true  line-of- 
sight  range-rate  between  a  GPS  satellite  and  receiver.  Deltarange  is  formed  in  the  receiver 
when  cycles  of  the  carrier-phase  are  counted  over  a  small  time  interval  divided  by  that 
time  interval.  This  time  interval  is  receiver  dependent,  typically  kept  small  compared 
to  the  measurement  output  rate.  For  the  Ashtech  P-XII  GPS  receiver,  deltarange  is 
calculated  over  a  1  millisecond  interval,  while  the  measurements  are  available  four  times 

^For  further  reference  on  GPS  measurements  and  applications  see  [10,19], 
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every  second  [19].  Since  carrier-phase  and  deltarange  measurements  are  observations  of 
the  same  signal,  they  are  highly  correlated. 

Occasionally,  one  sees  another  approximation  to  range-rate,  called  pseudo-deltarange 
or  pseudorange  deltarange.  Unlike  deltarange,  pseudo-deltarange  is  not  a  measurement 
available  from  a  GPS  receiver.  Instead,  pseudo-deltarange  is  formed  by  subtracting  succes¬ 
sive  pseudoranges  and  dividing  over  the  time  interval.  Time  intervals  of  pseudo-deltarange 
calculations  are  typically  on  the  order  of  1  second.  Because  pseudo-deltarange  uses  such  a 
long  time  interval  in  differentiating  position  information  into  velocity,  the  approximation 
to  range-rate  becomes  erroneous.  This  problem  can  be  exacerbated  in  the  case  of  high 
dynamics,  when  the  time  constant  of  the  dynamics  falls  close  to  or  below  the  Nyquist  rate. 
Since  the  deltarange  measurement  comes  from  the  receiver’s  carrier  loop,  it  is  independent 
of  the  pseudorange  measurement,  which  comes  from  the  code  loop.  Pseudo-deltarange, 
however,  is  clearly  correlated  with  pseudorange  measurements. 

1.3.2  Literature  Review.  Many  sources  of  information  were  drawn  upon  in  the 
preparation  of  this  research  including  numerous  papers  from  the  open  literature,  thesis 
reports,  textbooks,  course  notes  and  personal  interviews.  While  not  all  of  the  sources  are 
referenced  in  this  work,  the  referenced  sources  in  this  section  and  throughout  the  work 
provide  an  excellent  background  for  further  studies  in  this  area. 

As  of  this  writing,  the  only  paper  in  the  open  literature  describing  the  concept  of 
SARS  is  by  Raquet,  et  al  [36].  This  paper  also  presents  the  results  of  an  initial  test 
of  a  prototype  system  conducted  at  Holloman  AFB,  N.M.  Their  results  demonstrated 
the  feasibility  of  the  SARS  concept  in  determining  an  accurate  position  solution.  While 
SARS  is  unique  in  its  inverted  pseudolite  mode  of  operation,  there  are  other  trajectory 
determination  systems  that  use  GPS  methods  (carrier-phase  and  differential  C/A  code) 
for  kinematic  positioning.  These  are  described  by  Tang,  et  al,  [44],  Mahmood  and  Simpson 
[21],  Dougherty,  et  al  [8],  Robbins  [37]  and  Galijan  and  Gilkey  [9]. 

The  GPS  system  has  its  own  measure  of  velocity  (deltarange),  as  discussed  earlier. 
One  of  the  first  papers  to  quantify  GPS  velocity  accuracy  was  by  McGowan  [28].  May,  et 
al,  [23]  describe  the  properties  of  GPS  velocity  accuracy  as  well  as  applications  using  GPS 
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to  determine  velocities.  Methods  and  general  guidance  for  proper  treatment  of  deltarange 
measurement  in  a  positioning  system  are  given  in  Bach,  et  al  [1]  and  Brown  and  McBurney 
[3], 

The  previous  CIGTF-sponsored  thesis  work  of  Bohenek  [2],  Hansen  [11],  Mosle  [29], 
Negast  [35],  Solomon  [41],  Snodgrass  [40],  Stacey  [43]  and  Vasquez  [47],  provide  great 
insight  into  the  characterization  and  development  of  an  NRS,  as  well  as  the  background 
theory  on  inertial  navigation  systems,  error-state  extended  Kalman  filters,  GPS,  differential 
GPS  and  carrier-phase  GPS.  Also  of  interest,  Tobin  [46],  Chang  [5],  Schwartz  [42]  and 
Roecker  [38]  address  issues  concerning  adaptive  Kalman  filtering  and  target  tracking. 

1.3.3  Concept  of  Operation.  As  stated  previously,  the  design  of  SARS  departs 
radically  from  its  modern  predecessors,  CHAPS  and  CIRIS.  Unlike  these  reference  systems, 
SARS  does  not  use  an  INS.  Nor  will  SARS  require  aircraft  to  carry  a  large  payload  of 
sensors  and  measuring  devices.  Instead,  the  sensors  and  measuring  devices  wiU  be  located 
on  the  ground,  in  the  form  of  movable,  self-surveying  GPS  receivers.  Before  a  flight  test, 
engineers  will  be  able  to  lay  out  the  desired  locations  of  the  receivers  along  the  range, 
ensuring  adequate  receiver  geometry  and  coverage  during  the  flight. 

The  aircraft  will  be  outfitted  with  a  lightweight,  low  power  GPS  pseudolite,  designed 
to  be  mounted  on  aircraft  hard-points.  An  additional  GPS  pseudolite  will  broadcast  from 
a  known  position  on  a  mountain-top,  within  view  of  the  ground  receivers.  This  additional 
pseudolite  will  be  used  to  perform  a  between  receiver /pseudolite  double  difference  of  the 
carrier-phase  measurement.  This  double  differencing  method  eliminates  the  receiver  and 
pseudolite  clock  bias  as  well  as  reduces  errors  due  to  tropospheric  delay.  While  the  theory 
of  this  unique  carrier-phase  differencing  technique  is  still  under  development  [36],  this  thesis 
work  assumes  that  the  carrier-phase  measurements  will  yield  a  position  solution  with  an 
accuracy  of  0.1  m  la  3-D  RMS.  The  general  concept  of  operation  is  depicted  in  Figure  1.1. 

After  a  flight  test,  carrier-phase  position  data  will  be  collected  from  all  receivers. 
This  data  will  be  processed  to  resolve  the  receiver  clock  bias,  carrier-phase  ambiguities 
and  cycle  slips.  In  general,  a  cycle  slip  is  caused  by  a  loss  of  lock  of  the  carrier  signal 
between  the  GPS  transmitter  and  receiver.  Cycle  slips  may  be  caused  by  an  obstruction, 
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Figure  1.1  Concept  of  Operation 

intense  ionospheric  activity  or  receiver  power  problems  [19].  Since  SARS  is  limited  to  use 
in  the  troposphere,  no  cycle  slips  due  to  ionospheric  activity  are  anticipated. 

Cycle  slips  essentially  corrupt  the  future  time-history  of  the  carrier-phase  ambiguity 
term.  In  dynamic  applications,  cycle  slips  are  more  difficult  to  predict  due  to  the  additional 
Doppler  shift  created  by  the  aircraft’s  motion  [19].  Previous  AFIT  thesis  work  [2,11,29,47], 
has  addressed  the  problem  of  detecting  and  correcting  for  cycle  slips  in  highly  dynamic 
environments.  For  this  effort  it  is  assumed  that  the  carrier-phase  processing  will  correct 
the  position  solution  for  any  cycle  slips. 

1-4  Assumptions 

This  section  lists  the  assumptions  used  in  this  thesis  research.  These  assumptions 
are  defined  to  aid  the  reader  in  making  a  proper  evaluation  of  this  effort. 

1.  It  is  assumed  that  the  carrier-phase  data  has  been  processed  to  remove  ambiguities 
and  cycle  slips.  It  is  assumed  that  the  relevant  residual  carrier-phase  errors  include 
receiver  noise,  tropospheric  delay  and  multipath. 
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2.  It  is  assumed  that  the  processed  carrier-phase  solution  for  position  is  expected  to 
yield  errors  of  0.1  m  Icr  3-D  RMS.  This  position  error  is  assumed  to  include  the  error 
due  to  imprecise  knowledge  of  the  receiver  and  transmitter  (stationary  pseudolite) 
locations. 

3.  It  is  assumed  that  the  maximum  measurement  rate  of  the  GPS/pseudolite  receivers 
to  be  used  for  SARS  is  10  Hz. 

4.  It  is  assumed  differentially-corrected  C/A-code  pseudorange  and  deltarange  measure¬ 
ments  are  available. 

5.  It  is  assumed  that  all  GPS/pseudolite  receivers  are  perfectly  synchronized  in  time 
and  that  all  measurements  are  available  simultaneously. 

6.  Flight  trajectories  will  contain  a  range  of  dynamic  conditions  with  profiles  not  ex¬ 
ceeding  9  gees  and  10  m/s^  jerk.  These  profiles  are  assumed  to  be  representative  of 
the  dynamic  environment  of  modern  fighter  aircraft. 

7.  The  flight  profiles  used  in  this  thesis  effort  will  come  from  PROFGEN  [32].  MATLAB 
[27]  will  be  used  for  simulation,  data  reduction  and  visualization.  Explanations  of 
these  software  tools  can  be  found  in  Section  1.7.1. 

8.  Any  Monte  Carlo  analyses  to  be  conducted  wiU  be  the  results  of  50-run  simulations. 
Although  a  larger  number  of  Monte  Carlo  runs  would  produce  sample  statistics 
that  more  closely  reflect  the  true  underlying  error  statistics,  with  statistical  error 
asymptotically  approaching  zero  as  the  number  of  runs  approaches  infinity,  fifty  was 
decided  upon  to  keep  the  total  computation  time  within  reasonable  limits. 

1.5  Scope 

The  models  developed  in  this  effort  are  based  upon  information  obtained  from  the 
open  literature  and  from  sources  given  by  the  sponsor,  CIGTF.  The  velocity  algorithm 
target  specification  is  0.005  m/s  la  3-D  RMS,  at  a  measurement  rate  no  greater  than  10  Hz. 
Several  different  filtering  techniques  will  be  considered  in  simulation,  and  recommendations 
made  concerning  the  feasibility  of  a  real  world  implementation  of  these  techniques. 
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1.6  Approach/ Methodology 

Overall,  two  distinct  approaches  to  the  problem  will  be  investigated.  One  method 
involves  the  use  of  digital  filters  and  numerical  methods  techniques  to  evaluate  velocity 
as  the  derivative  of  discrete,  noisy  position  information.  The  second  technique  uses  the 
structure  of  the  Kalman  filter  to  estimate  velocity  given  a  simple  kinematic  model  and 
GPS  measurements  of  position  and  velocity. 

1.6.1  Digital  Filtering  with  Numeric  Differentiation.  The  digital  filtering  with 
numeric  differentiation  effort  is  a  straight-forward  approach  to  the  problem.  Given  noise 
corrupted,  discrete  position  measurements,  a  numerical  derivative  is  computed.  Two  tech¬ 
niques  for  implementing  a  numerical  derivative  are  explored.  First,  a  standard  central 
difference  approximation  to  the  derivative  is  derived  from  a  Taylor’s  series  expansion.  Sec¬ 
ond,  a  curve  fitting  algorithm  based  upon  cubic  splines  is  employed  to  fit  a  curve  to  the 
noise  corrupted  position  data.  An  analytical  derivative  of  the  fitted  curve  will  be  evaluated 
to  determine  velocity.  The  performance  of  the  two  techniques  will  be  compared.  Further¬ 
more,  digital  filtering  methods,  in  the  form  of  low  pass  filters,  are  used  for  noise  reduction 
both  before  and  after  the  derivative  operation. 

1.6.2  Kalman  Filtering.  The  Kalman  filtering  effort  will  devise  an  algorithm  for 
velocity  determination  for  two  different  measurement  strategies.  The  first  involves  using 
carrier-phase  data  that  has  been  processed  to  remove  ambiguities  and  the  effects  of  cycle 
slips.  The  second  strategy  involves  combining  unprocessed  pseudorange  and  deltarange 
measurements  to  determine  velocity.  Simulations  will  be  used  to  carry  out  this  research, 
and  simple  models  will  be  constructed. 

In  the  CHAPS  NRS,  an  error-state  extended  Kalman  filter  [2,11,39]  provides  esti¬ 
mates  of  the  reference  system  errors.  Figure  1.2  depicts  such  a  system.  The  estimated 
errors  are  subtracted  from  the  indicated  trajectory,  producing  the  estimated  reference  po¬ 
sition,  velocity  and  acceleration.  This  error-state  model  is  generated  from  higher  order 
truth  models  and  is  validated  for  performance. 
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In  order  to  create  a  proper  error-state  model  for  SARS,  a  truth  model  needs  to 
be  developed.  One  option  would  be  to  create  a  truth  model  based  upon  what  is  known 
about  carrier-phase  GPS  measurements  and  apply  it  to  this  “inverted”  GPS  problem. 
Such  truth  models,  however,  are  based  upon  a  mathematical  description  of  the  device,  and 
should  be  validated  by  empirical  data  obtained  from  the  “real  world.”  Since  SARS  is  still 
in  development,  it  is  unlikely  that  any  such  empirical  data  can  be  collected  at  this  point 
to  validate  such  a  model. 

It  is  possible  to  generate  an  error-state  model  based  upon  previous  carrier-phase 
models.  Such  a  model  would  need  to  be  modified  based  upon  our  assumptions  of  the 
operating  environment.  After  the  model  is  modified,  it  should  be  validated.  This  validation 
would  require  significant  resources  [39]  beyond  the  scope  of  this  thesis  effort. 

Even  without  formal  validation  it  is  possible  to  use  a  model;  for  example,  one  mod¬ 
ification  might  involve  removing  model  states  corresponding  to  ionospheric  errors,  since 
the  GPS  signal  will  not  be  traveling  through  the  ionosphere.  In  this  case,  it  appears  that 
our  assumptions  about  the  operating  environment  will  improve  the  accuracy  of  our  results. 
On  the  other-hand,  new  sources  of  error  may  exist  and  should  be  introduced  to  the  model. 
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Multipath  error,  while  not  considered  in  past  AFIT  carrier-phase  research  [2,11],  may  be 
a  significant  source  of  error  for  SARS.  If  so,  what  effect  does  altitude  have  on  multipath 
error?  Does  the  multipath  portion  of  the  model  need  to  be  adjusted  for  altitude  changes, 
and  if  so,  how?  While  it  may  be  prudent  to  consider  the  modeling  of  multipath  as  a  func¬ 
tion  of  altitude,  it  is  difficult  to  predict  the  impact  an  altitude  change  should  have.  Hence, 
without  real  world  validation,  model  confidence  is  quickly  lost. 

A  heuristic  approach  to  the  model  formulation  is  suggested.  If  one  assumes  that 
a  valid  error-state  model  cannot  be  confidently  constructed,  then  the  only  information 
remaining  to  be  exploited,  aside  from  the  carrier-phase  data  itself,  is  the  fact  that  the 
aircraft  will  be  in  controlled  flight,  along  some  trajectory.  If  an  aircraft  trajectory  can  be 
broken  down  into  distinct  phases,  the  kinematics  can  be  represented  by  simple  models, 
such  as: 

•  Constant  Velocity,  Acceleration  modeled  as  a  white  Gaussian  noise 

•  Constant  Acceleration,  Jerk  modeled  as  a  white  Gaussian  noise 

•  Constant  Acceleration,  Jerk  modeled  as  a  1st  order  Markov  process 

It  is  possible  that  one  of  these  models  may  be  generic  enough  to  provide  an  adequate 
velocity  estimate  over  the  entire  flight,  or,  each  model  may  be  valid  only  during  a  particular 
segment  of  a  flight.  In  this  case,  an  adaptive  filter  could  be  devised  that  changes  in  response 
to  flight  conditions. 

1.7  Materials,  Data  and  Equipment 

No  special  provisions  are  required  for  this  research.  The  Sun  SPARC  network  in  the 
AFIT  Guidance  and  Navigation  Laboratory  as  well  as  the  XR5M  6-channel  and  12-channel 
GPS  receivers  will  be  used  for  experimentation  during  this  thesis. 

1.7.1  Simulation  Software. 

1.7. 1.1  PROFGEN.  The  Avionics  Directorate  of  Wright  Laboratory  devel¬ 
ops  and  maintains  a  software  package  called  PROFGEN  (PROFile  GENerator).  PROF- 
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GEN  allows  for  the  creation  of  flight  trajectories  of  a  vehicle  under  continuous  control  over 
the  Earth’s  surface.  This  is  accomplished  by  solving  the  equations  of  motion  of  a  zero- 
mass  body  responding  to  maneuver  commands  specified  by  the  user.  These  maneuvers 
can  consist  of  vertical  turns,  horizontal  turns,  rolls,  straight  flights  and  sinusoidal  heading 
changes.  More  detailed  information  is  available  in  the  PROFGEN  documentation  [32]. 

1.7. 1.2  MATLAB.  MATLAB  [27]  is  a  technical  computing  environment 
for  high  performance  numeric  computation  and  visualization.  MATLAB,  which  stands  for 
matrix  laboratory,  is  useful  for  conducting  simulation  and  analysis  of  linear  (and  non-linear) 
systems.  MATLAB  features  a  powerful  matrix  interpreter,  a  user-friendly  environment, 
a  variety  of  reliable  numerical  integration  algorithms  and  a  graphical  simulation  environ¬ 
ment  called  SIMULINK.  MATLAB,  SIMULINK  and  several  “toolboxes”  were  used  in  the 
construction  of  the  simulations  and  their  analysis. 

1.8  Overview  of  Thesis 

The  methodology  in  this  chapter  is  meant  to  motivate  the  development  of  algorithms 
to  determine  the  velocity  of  an  aircraft  flying  in  a  trajectory.  CIGTIF’s  performance  objec¬ 
tive  for  SARS  is  a  velocity  accuracy  of  0.005  m/s  3-D  RMS  Icr  at  a  maximum  measurement 
rate  of  lOHz.  Two  approaches  will  be  considered:  a  digital  filtering  with  numeric  differen¬ 
tiation  approach  and  a  kinematic  model-based  Kalman  filter  approach. 

This  thesis  consists  of  five  chapters.  Chapter  I  has  presented  the  problem  to  be  solved 
and  the  approach  to  be  taken  in  solving  it.  Chapter  II  covers  the  theory  involved  in  this 
thesis.  Topics  covered  in  Chapter  II  include  an  overview  of  numerical  differentiation,  GPS 
measurement  equations,  digital  filtering  and  Kalman  filter  theory.  The  theory  involved 
with  GPS  measurements  includes  how  the  measurement  equations  are  derived,  in  light  of 
the  inverted  pseudolite  concept.  Chapter  III  presents  the  noise  models  and  Kalman  filter 
models  implemented  in  this  thesis.  The  results  of  the  data  evaluations  and  simulations  are 
presented  in  Chapter  IV  and  Chapter  V  gives  conclusions  and  recommendations  for  future 
research. 
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II.  Theory 


2.1  Overview 

This  chapter  describes  the  general  theory  pertaining  to  this  research.  The  first  sec¬ 
tion  outlines  the  implications  of  noise  corrupted  discrete-time  position  measurements  in 
velocity  determination.  The  GPS  measurement  section  presents  the  equations  involved 
with  the  three  basic  GPS  measurements:  pseudorange,  carrier-phase  and  deltarange;  fur¬ 
ther  information  on  GPS  theory  is  available  in  [10].  The  Kalman  filtering  section  provides 
an  overview  of  Kalman  filtering/smoothing.  Readers  unfamiliar  with  these  Kalman  filter¬ 
ing  topics  are  referred  to  [24,25].  Finally,  the  last  section  presents  the  types  of  digital 
filters  used  in  this  research. 


2.2  Approximating  Velocity 


GPS  measurements  of  the  aircraft  position  will  be  available  at  discrete  time  intervals. 
These  measurements  can  be  solved  in  Earth  Centered  Earth  Fixed  (ECEF)  coordinates 
and  expressed  as  follows: 

[  xiti)  1 


r{ti)  = 


y{ti) 


(2.1) 


z{ti) 


The  velocity  of  the  aircraft  can  then  be  defined  as  the  derivative  of  position: 


r(t.)  s  lim 


r(ti  +  h)-  r(ti) 
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(2.2) 


Unfortunately,  this  equation  cannot  be  properly  evaluated  as  h  is  limited  to  the  sampling 
rate  of  SARS. 


2.2.1  Taylor’s  Series  Approximation  to  the  Derivative.  One  method  used  to  ob¬ 
tain  the  numerical  derivative  approximation  and  maintain  stability  is  a  truncated  Taylor’s 
series  expansion.  The  Taylor’s  series  for  a  function  r{x)  at  {x  +  h)  expanded  about  x  is: 

/  ^  //  M  r"(x)h?  T"'(x)h^  . 

r{x  +  h)  =  r(x)  +  r  (x)h -\ - — — -| - - \-h.o.t.  (2.3) 
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Similarly  for  r{x)  at  {x  —  h)  we  have 


T"(x)h?'  r'"(x)h^ 

t{x  —  h)  =  r{x)  -  r  (x)h -\ - j - j - \- h.o.t.  (2-4) 

Combining  (2.3)  and  (2.4)  we  have: 

[r(.  +  k)  -  r(x  -  ft)]  =  2rXx)  +  2’^^  +  2^^  +  k.o.t.  (2.5) 


Notice  that  the  even  numbered  derivative  terms  are  eliminated.  Truncating  (2.5)  to  first 
order  and  rearranging  yields  the  classic  central  difference  approximation  of  the  first  deriva¬ 
tive: 
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2.2.2  Extrapolating  the  Central  Difference  Approximation.  In  (2.5),  the  even 
order  derivatives  in  the  Taylor’s  series  cancel  out  fortuitously.  It  is  possible  to  combine 
additional  equations  to  cancel  out  additional  higher  order  terms.  Consider  the  following 
equation: 

[r{x  +  2h)  -  r(x  -  2h)]  =  4r'{x)  +  2  •  2^-M-  +  2  •  +  h.o.t.  (2.7) 


which  can  be  combined  with  a  scaled  version  of  (2.5)  to  remove  the  third  order  derivative: 

( x')h^ 

2"[r(a;  +  h)-  r{x  -  h)]  -  [r(a;  +  2h)  -  t{x  -  2h)]  =  12r'(a:)  -f  2  •  (2'*  -  2®)—^^  +  h.o.t.  (2.8) 


Neglecting  the  fifth  order  and  higher  derivatives,  we  can  rearrange  (2.8)  to  solve  for  the 
following  central  difference  equation: 


r'(a;)  ^  8 


r(a;  +  h)  —  t{x  —  h)  r{x  -f  2h)  —  r{x  —  2h) 


I2h 


12h 


(2.9) 


This  extrapolation  method  can  be  pursued  further  to  generate  approximations  to  deriva¬ 
tives  with  better  residual  error  characteristics,  at  the  expense  of  requiring  additional  data 
points.  The  central  difference  equation  approximation  for  the  first  derivative  is  well  docu¬ 
mented  [6,15,22,45]  in  numerical  methods  texts.  Equations  (2.10)  ,(2.11)  and,  (2.12)  show 
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the  algorithms  tested  in  this  research. 


r'(^) 

r'(.T) 

r'(a:) 


^  t{x  +  h)-  r{x  -  h) 

2h 

^  ^  r{x  +  h)  -  r{x  -  h)  r{x  +  2h)  -  r(a;  -  2h) 

m  m 

r(a;  +  h)  -  r(a:  -  h)  ^  r(x  +  2h)  -  r{x  -  2h) 

60h  60h 

r{x  +  3h)  -  r(x  —  3/i) 

60h 


(2.10) 

(2.11) 


(2.12) 


For  best  precision,  the  numerical  differentiation  algorithm  would  be  extrapolated  as  far  as 
higher  order  derivatives  exist  for  the  function  in  question.  Due  to  the  sampled  data  nature 
of  the  system  under  consideration  in  this  research,  there  is  no  guarantee  that  the  function 
is  continuous,  especially  in  areas  of  high  dynamics. 


2.2.3  Error  Analysis  of  the  Velocity  Approximation.  There  are  many  contributing 
factors  in  the  overall  error  of  the  velocity  approximation  algorithms.  Three  sources  of  error 
in  the  velocity  approximation  are  considered:  error  due  to  digital  sampling,  error  due  to 
truncation  of  the  Taylor’s  series  and  error  due  to  corrupted  position  data.  This  section  is 
concerned  with  evaluating  the  extent  of  these  errors. 


2.2.3. 1  Error  Due  to  Digital  Sampling.  If  we  treat  the  position  measure¬ 
ments  as  digital  sampling  of  a  continuous  time  signal,  one  must  be  concerned  with  the 
phenomena  of  aliasing.  In  the  frequency  domain,  aliasing  results  in  higher  frequencies  of  a 
signal  impersonating  lower  frequencies.  In  the  time  domain,  aliasing  results  in  corruption 
of  the  discrete-time  signal.  In  order  to  avoid  aliasing,  it  is  necessary  to  sample  at  a  rate 
at  least  two  times  higher  than  the  highest  frequency  of  the  input  signal.  This  is  known  as 
the  Nyquist  rate. 

To  determine  the  highest  frequency  component  that  the  input  signal  may  contain, 
we  examine  the  case  of  a  high  performance  fighter  aircraft  in  a  demanding  flight  profile. 
In  Section  4.3.1,  analysis  confirms  that  Nyquist  criterion  is  met  for  the  sampling  rates  of 
interest  in  this  research. 
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2.2.3. 2  Error  Due  to  Truncation.  Neglecting  higher  order  terms  in  the 
Taylor’s  series  will  cause  errors  of  magnitude  h"  where  h  is  the  sampling  interval  and  n  is 
the  order  of  the  derivative  of  the  first  neglected  term  [6,15,22].  For  example,  the  neglected 
terms  in  (2.10)  are: 


2  •2-’ 


3! 


+  2-2'’ 


5! 


+  h.o.t. 


(2.13) 


and  the  order  of  the  error  is  h^.  Similarly,  (2.11)  and  (2.12)  are  of  order  and  hJ 
respectively. 


2. 2.3. 3  Implications  of  Noisy  Position  Data.  Now  consider  errors  in  the 
position  information.  Each  coordinate  is  corrupted  by  some  error  S: 


r(i.-) 


x{ti)  +  6x{ti) 

y{ti)  +  hiu) 

z{ti)  +  6z{ti) 


(2.14) 


where  x{ti)^  y(ti)  and  z{ti)  represent  true  position.  Using  (2.6)  the  velocity  equation  for 
the  X  direction  becomes: 


dx{ti)  ^  [a:(t,+i)  +  ^a:(t,+i)]  -  [a;(ti_i)  +  ^a:(t.-i)] 

dt  ~  2At 

xjtj+i)  -  x{ti_-i)  6x(ti+i)  -  6x{ti^i) 

2At  2At 


(2.15) 


The  velocity  equation  can  be  expressed  as  a  sum  of  two  terms:  the  computed  double  sided 
velocity  approximation  and  the  error  in  computed  velocity.  In  order  to  generate  accurate 
velocity  estimates,  it  is  important  to  minimize  the  effect  of  the  error  term  in  Equation 
2.15.  One  way  to  make  this  term  approach  zero  is  to  make  the  sampling  period  very  large 
compared  to  the  magnitude  of  the  difference  in  the  error.  This  solution  is  not  practical 
since  large  sampling  periods  invalidates  the  difference  equation  as  an  approximation  to 
velocity. 

Another  way  to  reduce  the  effect  of  the  error  term  is  to  make  the  numerator  much 
smaller  than  the  denominator.  If  the  errors  ^a:(t,+i)  and  5a;(t,_i)  are  time  correlated,  they 
may  cancel  each  other  out,  depending  on  the  sampling  rate  and  the  time  constant  of  the 
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errors.  Unfortunately  there  are  practical  limits  to  the  sampling  rate.  If  there  is  a  white 
component  to  the  noise,  or  if  the  error  terms  are  not  exactly  correlated,  there  will  be  a 
non-zero  remainder.  As  the  sampling  rate  is  increased,  this  remainder  becomes  amplified 
by  the  denominator  term. 

Care  must  be  taken  when  choosing  a  sampling  rate;  a  trade-off  between  the  accuracy 
of  the  velocity  approximation  and  the  amount  of  noise  amplification  must  be  made.  In 
many  applications,  such  as  navigation,  it  is  prudent  to  eliminate  as  much  noise  as  possible 
from  the  position  measurements  before  using  the  differencing  equation  to  approximate 
velocity. 

2.3  Global  Positioning  System  Measurements 

2.3.1  Pseudorange.  The  basic  pseudorange  observable  is  a  difference  between 
the  time  of  transmission  of  the  GPS  signal  and  the  time  of  arrival,  based  upon  a  code 
correlation  technique.  The  pseudorange  observation  equation  given  by  [10]  is: 


P  — II  Rj  II  "fC  •  bfc)  -|-  bfrop  -|-  bjnpath  T  ^ion 


where:  p 

= 

is  the  measured  pseudorange 

II  ~  II 

the  true,  but  unknown,  range 

Ri 

ECEF  coordinates  of  the  i*’'  receiver 

Rj 

ECEF  coordinates  of  the  pseudolite 

6pc 

the  pseudolite  clock  offset  from  true  GPS  time 

^rc 

the  receiver  clock  offset  from  true  GPS  time 

C 

speed  of  light 

^ion 

= 

the  range  equivalent  ionospheric  delay  term 

^trop 

= 

the  range  equivalent  tropospheric  delay  term 

^mpath 

= 

the  error  due  to  multipath  effects 

Since  the  pseudolites  used  for  this  application  will  not  transmit  signals  through  the 

ionosphere,  the  error  term  Sion  may  be  eliminated.  The  resulting  pseudorange  equation  is: 

p  — II  Rj  Ri  |j  -j-C  ■  (^pc  ^rc)  "t"  ^trop  T  ^mpath  (2.16) 
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which  is  valid  for  both  the  airborne  and  ground-based  pseudolites.  The  receiver  location, 
Ri,  is  known  a  priori  since  the  ground  receivers  use  carrier  phase  surveying  techniques  to 
determine  precisely  their  position.  Similarly,  the  ground-based  pseudolite  position  is  also 
known  a  priori.  Thus,  a  good  estimate  of  the  range  between  the  ground-based  pseudolite 
and  receiver  is  known.  The  true  range  can  be  expressed  as  a  sum  of  two  terms,  the  surveyed 
or  estimated  range,  ||  Rj  —  Ri  ||,  and  the  error  in  the  survey,  Sr. 

11  Rj  ~  Ri  11=11  Rj  —  Ri  II  (2-17) 

By  subtracting  the  known  quantity  1|  Rj  —  Ri  ||  from  both  sides  of  the  ground-based 
pseudorange  observation  (2.16),  we  define  a  differential  correction  term  Ap: 

Ap  =  Pgnd-  \\  Rj  -  Ri  \\ 

—  [c  '  (^pc  ^rc)  H"  ^trop  (2.18) 

where  the  subscript  gnd  denotes  terms  corresponding  to  the  ground-based  pseudolite. 

The  quantity  expressed  in  (2.18)  can  be  subtracted  from  the  pseudorange  observation 
equation  for  the  airborne  pseudolite,  providing  a  differentially  corrected  pseudorange,  p*. 

P*  -  Pair  "  Ap 

—  II  Rj  Ri^  II  T  [c  '  (^pc  ^rc)  ^trop  “I”  ^mpath^Q^r 

[c  *  (^pc  ^rc)  H“  ^trop  “b  ^rnpath\gyif^ 

~  II  Rj  Ri  1!  “b  '  ^pc  H“  ^trop  H"  ^mpaihlfiiy, 

[c  *  Sp(^  Sfj-Qp  +  ^mpathlgnd  (2.19) 

where  the  subscript  air  denotes  terms  corresponding  to  the  airborne  pseudolite.  This 
eliminates  the  common  receiver  clock  bias  since  [^rc]^„d  “  \^rc\air->  reduce  the 

tropospheric  delay  effect. 

To  further  reduce  the  error  in  (2.19),  a  tropospheric  error  correction  can  be  calcu¬ 
lated.  In  GPS  applications,  the  tropospheric  delay  is  considered  a  function  of  humidity. 
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temperature  and  altitude  at  the  receiver.  The  implications  of  modeling  the  tropospheric 
delay  term  is  developed  in  Section  3.3.1. 

2.3.2  Carrier-Phase.  This  section  develops  the  essential  pseudolite-based  carrier- 
phase  equations  required  for  this  effort.  For  the  reader  unfamiliar  with  carrier-phase  meth¬ 
ods,  a  more  general  derivation  of  the  GPS  carrier-phase  equations  is  presented  in  [2, 11]. 

A  carrier-phase  measurement  is  derived  from  the  phase  shift  between  the  transmitter¬ 
generated  signal  and  receiver-generated  signal.  In  GPS  applications,  the  transmitter  is  a 
satellite;  in  this  application,  the  transmitter  is  one  of  two  pseudolites.  One  pseudolite  is 
located  at  a  surveyed  location  on  the  Earth’s  surface,  in  view  of  all  the  receivers.  The 
second  pseudolite  is  mounted  on  an  aircraft,  flying  in  a  trajectory  over  the  range. 

The  carrier-phase  observation  equation  [2]  for  GPS  pseudolites  is  given  by: 

$  =  — /  •  (^pc  —  Kc)  —  ~  •  (II  Rj  ~  R-i  II  ~^ion  +  ^trop  +  ^mpath)  (2.20) 

where:  $  =  the  phase-range  in  cycles  between  pseudolite  and  receiver 

/  =  the  frequency  of  transmission 

Spc  =  the  pseudolite  clock  offset  from  true  GPS  time 

Src  =  the  user  clock  offset  from  true  GPS  time 

c  =  speed  of  light 

II  Rj  -  Ri  II  =  the  true,  but  unknown,  range 

Sion  —  the  range  equivalent  ionospheric  delay  term 

Strop  —  the  range  equivalent  tropospheric  delay  term 

^mpath  =  the  error  due  to  multipath  effects 

Since  the  pseudolite  signals  do  not  pass  through  the  ionosphere,  the  term  Sio„  can  be 
eliminated  from  (2.20).  The  resulting  carrier-phase  observation  equation  for  pseudolites  is: 

^  =  —f  •  {hpc  —  Src)  —  ^  ■  (II  Rj  ~  R>  II  +*^trop  +  ^mpath)  (2.21) 

Because  the  carrier-phase  measurement  is  a  phase  shift,  it  represents  only  a  fraction 
of  a  wavelength  of  the  carrier.  The  phase-range  measurement  at  some  time  epoch,  t,  is 


represented  by  the  following  equation: 


$(t)  =  ^irac{t)  +  +  N  {U)  (2.22) 

where  $/rac(^)  is  the  fractional  part  of  the  total  wavelength,  t)  is  an  integer  number 

of  phase  cycles  from  an  initial  epoch,  to  the  current  epoch,  t,  and  N{to)  is  an  integer 
phase  ambiguity  term.  The  phase  ambiguity  term  is  also  known  as  the  cycle  ambiguity 
and  it  represents  the  difference  between  the  true  integer  count  at  time  to,  and  the  current 
integer  count  at  tg  measured  or  calculated  by  the  receiver  [11]. 

The  carrier  phase  measurement,  ^measured{t)i  is  equal  to  the  sum  of  the  fraction 
observation  at  time  epoch,  t,  and  the  integer  count  at  the  same  time  epoch,  t,  and  can  be 
represented  by: 

^measuredil)  —  ^/rac(f)  T  (2.23) 

The  total  phase  range  at  time  epoch,  t,  from  (2.22)  can  now  be  written  as: 

$(i)  =  $m.a.«re<i(<)  +  iV(t)  (2.24) 

Substituting  (2.21)  into  (2.24)  produces  the  measured  range  for  the  pseudolite  carrier-phase 
observable: 

^measured(t)  —  ^(0  -^(0 

—  ~f  '  (^pc  ~  ^rc)  —  ~  ■  (II  II  -i-^trop  +  bypath)  —  N (t)  (2.25) 

Multiplying  (2.25)  by  the  negative  of  the  cycle  wavelength,  -X,  converts  from  cycle  units 
to  distance  units: 

^(t)  —  X  •  $measure<i(t) 

=  II  Rj  —  Rj  II  -f-C  •  {6po  —  Src)  +  btrop  +  bmpath  +  X  ■  N (t)  (2.26) 

2.3.2. 1  Carrier  Phase  Double  Differencing.  The  between-receiver/satellites 
double  differencing  method  described  in  [2,36]  is  the  method  to  be  implemented  in  this 
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research.  It  is  modified  slightly  to  a  between-receiver/between-pseudolites  method  and  is 
depicted  in  Figure  2.1.  This  method  subtracts  a  between-receiver  single  difference  with 


i-th  Ground  Based  Receiver  (GBR)  j»th  Ground  Based  Receiver  (GBR) 


Figure  2.1  Illustration  of  Between  Receivers/Pseudolites  Double  Difference 

another  between-receiver  single  difference  using  the  same  receivers  and  the  ground-based 
reference  pseudolite.  Additionally,  the  integer  ambiguity  term  must  be  estimated.  In  [36], 
preliminary  results  using  this  double  differencing  method  with  conventional  integer  ambi¬ 
guity  resolution  techniques  are  presented.  The  result  of  this  double  differencing  method 
is  the  elimination  of  pseudolite  and  receiver  clock  error  biases.  An  atmospheric  error 
term  reduction  is  dependent  upon  the  relative  magnitudes  and  directions  of  the  ranges 
between  the  transmitters  and  receivers  as  well  as  the  homogeneity  of  the  atmosphere,  and 
subsequently,  is  highly  application  specific. 

2.3.3  Deltarange.  By  differentiating  the  carrier-phase  measurement  equation 
(2.26),  we  can  derive  an  approximation  to  the  range  rate: 

^  f  A  II  Rj  —  Rj  II  +C  ■  A{6pc  —  6rc)  +  ^^trop  +  ^mpath  +  ^  '  AA (t) 
dt  At->-o  \  At 

which  approaches  the  true  range  rate  as  At  approaches  zero.  The  true  range  rate  equation 
is  given  by: 

d  II  Rj  Rj  II  dSjj.Qp  dbpf^  dSfc 

It  =  — It —  + 
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d  II  Rj  —  Ri  II  dS 

error 

dt  dt 


(2.28) 


where  terror  =  ^trop  +  ^mpath  +  ^pc  ”  ^rc-  This  assumes  that  the  integer  phase  ambiguity 
term,  N{t),  remains  constant  over  the  time  interval  such  that  AN{t)  =  0. 

2. 3. 3.1  Solving  the  Range  Rate  Equation  for  User  Velocity.  User  velocity 
is  embedded  in  the  first  term  of  (2.28).  Evaluating  this  term  we  have: 


d  II  Rj  —  Ri 
dt 


+  (%  -  ViY  +  -  ZiT 

{Xj  -  Xi){xj  -  Xj)  +  {yj  -  yi){yj  -  y,)  +  {zj  -  Zi){zj  -  i.) 
\/{xj  -  XiY  +  (yj  -  yif  +  (Zj  -  Zi^ 


(2.29) 


where  i  denotes  the  i**  pseudolite  and  and  j,  the  j*'*  receiver.  Letting  Oi,  a^,  and  03  be 
the  directional  derivatives: 


Oi  = 


as  = 


Oia  = 


/ixj-Xi)^+(yj-yiy+(zj-Ziy 

iy,-yi) 


/(xj-Xi)^+(yi-yi)^+(zj-Ziy 


Equation  (2.29)  can  be  expressed  as: 


d  II  II  ,  .  .  .  ,  .  . ,  ,  .  . , 

- -  Xi)  +  aaiyj  -  yO  +  asizj  -  Zi) 


(2.30) 


Substituting  (2.30)  into  (2.28)  yields: 


—  ==  ai{xj  -  Xi)  +  aaiyj  -  Vi)  +  03(2^  -  Zi)  + 


Rearranging  this  equation  to  separate  pseudolite  velocity  from  user  velocity: 


—  +  aiXi  +  aaVi  +  a3ii  =  a^Xj  +  aaVj  +  a^Zj  +  c — 


(2.31) 


(2.32) 
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For  the  roving  pseudolite  on-board  the  aircraft,  the  velocities  of  the  receivers  are  zero 
with  respect  to  the  ECEF  frame.  Equation  (2.32)  for  the  airborne  pseudolite  becomes: 

.  .  .  dSerror  /„  „„n. 

—  =  aiXj  -1-  a-iVj  -h  a^Zj  -f  c  -  -  (2.33) 

dt  air  di  air 

For  the  ground-based  pseudolite,  not  only  are  the  velocities  of  the  receivers  zero,  but  also 
the  velocity  of  the  pseudolite  itself  is  zero.  Equation  (2.32)  for  the  ground-based  pseudolite 
becomes: 

dt  grid 

If  we  assume  that  simultaneous  measurements  of  the  pseudolites  is  possible,  the 
ground-based  pseudolite  can  be  used  as  a  differential  correction  to  the  airborne  pseudolite. 
Since  the  same  receiver  is  used  to  detect  both  pseudolite  signals,  the  quantities 
and  are  equal.  For  the  remaining  error  terms,  c^^^andc^^^^^,  it  is  not 

known  what  benefit  differential  corrections  will  have.  Assuming  the  pseudolites  are  of 
the  same  make  and  model,  they  should  use  the  same  clock  device  and  have  similar  drift 
characteristics,  However,  the  airborne  clock  will  be  subjected  to  higher  g-forces  and 
may  suffer  from  the  effects  of  g-induced  disturbances  to  the  oscillator  [20]. 

Nonetheless,  the  tropospheric  error  rate,  may  be  small  enough  over  such  a 

short  period  of  time  to  neglect  entirely.  This  assumes  that  there  is  a  high  correlation  in 
tropospheric  error  between  subsequent  measurements,  which  is  likely  given  that  sequential 
signals  are  passing  through  (relatively)  the  same  atmosphere.  While  the  pseudolite  signals 
are  travelling  through  essentially  the  same  atmosphere  (wet  and  dry  refract! vities),  the 
tropospheric  error  terms  are  also  dependent  on  user  altitude,  range  and  elevation  angle 
[17].  Experimentation  may  be  useful  in  determining  the  usefulness  of  such  a  differential 
correction.  The  differentially  corrected  range  rate  equation  is: 
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dt  diff 


dt  air 


dt  gnd 


.  ,  .  ,  .  ,  ^^error  d^error 

=  aiXj  +  a2j/j  +  asZj  +  c  -  c— — - 

Cut  air  (Vt  gnd 


=  aiXj  +  Q!2?/j  +  aaZj  +  c- 


dS, 


trop 


+  C 


d6. 


mpaih 


dt  air— gnd  dt  air  — gnd 


+  c~ 


dS, 


pc 


—  diXj  +  0:2?/j  + 


db^ 


dt  air— gnd 


dt  air  — gnd 

(2.35) 


2.4  Kalman  Filtering/Smoothing 

Although  solutions  to  the  GPS  ranging  equations  are  non-linear,  in  this  exploratory 
research,  the  Kalman  filter  measurement  models  were  kept  linear.  This  is  accomplished  by 
solving  the  ranging  equations  before  presenting  them  to  the  Kalman  filter.  This  was  done 
to  keep  the  analysis  focused  on  assessing  the  performance  and  suitability  of  the  Kalman 
filter  algorithm  itself.  If  the  desired  level  of  performance  can  be  achieved  with  the  standard 
Kalman  filter,  and  possibly  exceeded,  the  study  of  the  effects  of  compensation  for  non-linear 
measurement  equations  is  warranted. 

While  the  Kalman  filter  lends  itself  to  use  in  real-time  applications,  it  is  also  useful  in 
non-real  time  estimation  problems.  In  fact,  by  taking  advantage  of  information  contained 
in  “future”  measurements,  a  better  state  estimate  can  be  provided  in  most  cases  [25].  This 
use  of  “future”  measurements  in  a  Kalman  filter  structure  describes  the  optimal  smoothing 
problem. 

While  many  different  smoother  formulations  exist,  the  fixed-interval  smoother  as  de¬ 
scribed  in  [25],  is  most  applicable  to  this  research.  The  fixed-interval  smoother  is  especially 
useful  for  post-mission  data  analysis  as  it  allows  all  measurement  data  collected  to  be  used 
to  generate  the  state  estimates  at  a  given  time,  f,-. 

Considering  a  discrete-time  model  of  states  x(f,),  and  measurements  z{ti): 

x(ti+i)  =  ^(,ti+i,ti)x{ti)  +  Bd{ti)u{ti)  +  (2.36) 

z{ti)  =  H(ti)x(ti)  -t-  \i{ti)  (2.37) 
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where  the  initial  conditions  x(to)  have  a  mean  Xq  and  covariance  Pq.  States  are  propagated 
by  a  discrete-time  state  transition  matrix  Also,  Wd(t,)  is  a  discrete-time  white 

Gaussian  noise  with  covariance: 


Qd(<*)  ti  =  tj 

0 


(2.38) 


and  is  independent  of  the  measurement  noise  vector,  v(t,),  which  also  is  described  as  zero- 
mean  white  Gaussian,  with  noise  covariance,  R(tj): 


R(t.)  ti  -  tj 

0 


(2.39) 


The  smoothing  procedure  can  be  thought  of  as  having  three  steps.  First  a  con¬ 
ventional  forward  running  Kalman  filter  generates  state  estimates,  denoted  x(tf ),  where 
the  (+)  superscript  indicates  the  state  estimate  after  a  measurement  update  at  time  t,. 
Second,  a  backward  running  Kalman  filter,  independent  of  the  forward  running  filter,  com¬ 
putes  state  estimates  Xj(t~),  where  the  (-)  superscript  indicates  the  state  estimate  before 
a  measurement  update  at  time  ti.  Third,  the  state  estimates  are  combined  optimally  to 
produce  :k{tiftf),  by  considering  them  as  two  separate  observations  of  x(tj),  weighting  each 
by  factors  related  to  the  covariances  P{tf)  and  Pj,(t,“).  Computationally,  it  is  more  effi¬ 
cient  to  run  and  store  the  results  of  the  backward  filter  first,  then  combine  steps  one  and 
three  into  a  single  algorithm  that  computes  the  forward  filter  estimate  and  the  smoothed 
estimate  simultaneously. 

The  fixed-interval  smoothing  algorithm,  as  presented  in  [25],  is  shown  below.  A 
conventional  forward  running  Kalman  filter  with  propagation  equations: 

M^k+i)  =  ^{tk+uh)^t+)  +  Bd{ti.)u{tk)  (2.40) 

P(4-n)  =  ^(4  +  l?4)P(tfc  )^"^(it-|-l5tj:)  +  Gr(i(4)Qd(tjS;)Grd^(/)i.)  (2-4l) 
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and  measurement  update  equations: 


K(4)  =  P(4-)HT(4)  [H(4)P(4“)H'^(4)  +  R(4)]"'  (2.42) 

x(4)  =  x(4)  +  K(4)  [zfc -H(4)x(4)]  (2.43) 

P(t+)  =  P(4-)-K(4)H(4)P(4)  (2.44) 

and  initial  conditions: 

x(4)  =  xo  (2.45) 

P(4)  =  Po  (2.46) 


are  used  to  generate  a  state  estimate  and  covariance  for  each  time.  A  backwards  run¬ 
ning  Kalman  filter  is  implemented  in  inverse-covariance  form  since  no  a  priori  statistical 
information  is  known  about  the  initial  conditions.  Measurement  updates  are  given  by: 

y6(4+)  =  yi(4-)  +  H'^(4)R-'(4)z.  (2.47) 

Pr'(4+)  =  Pr'(4-)  +  HT(4)R-i(4)H(4)  (2.48) 

The  states  are  propagated  backwards  in  time  by: 

J(4)  =  Pr'(4+)Gd(4-i)  [G/(4-i)Pr'(4+)Gd(4-i)  +  Qr'(4-i)] 


L(4)  -  I- J(4)g7(4-.i) 

yi(4_J  =  -^^(4,4+i)L(4)[y6(4+)-Pr'(4'')Bd(4-i)u(4-i)]  (2.49) 

PrH4"-i)  =  $"'(4,4-n){L(4)Pr'(4+)L'^(4) 

+  J(4)Q(i  (2.50) 

with  initial  conditions: 

yi(tf)  =  0  (2.51) 

Pr'(4“)  =  0  (2.52) 
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The  smoothed  estimate  is  generated  by  combining  x(t+),  P(t+),  yi{ti  ),  and  Pj  in: 

XiU)  =  [l  +  P(t+}Pr^(tr)]“^ 

W{U)  =  P(t+)XT(i,) 

Y{t,)  =  i-w(t,)prHir)w^O.-) 

P(t,7t;)  =  Y{U)P{tt)Y'^{ti)  +  W{ti)Pi-\t-)W^{U)  (2.53) 

±{ti/tf)  =  X{ti)x{tf)  +  P{ti/tj)yi(tr)  (2.54) 

where  X{ti),  W(t,)5  a-nd  Y(t,)  are  defined  as  intermediate  terms  used  to  evaluate  the 
smoothed  state  estimate,  and  associated  covariance,  P{tiftf).  The  notation 

{tiftf)  stands  for  “at  time,  tj,  based  on  measurements  through  time  t,  and  on  measurements 
from  time  tj". 

A  metric  has  been  developed  to  determine  the  relative  benefit  of  Kalman  smoothing 
over  Kalman  filtering  alone.  This  metric,  the  Fraser  smoothability  criterion  [25],  is  based 
upon  an  eigenvalue  decomposition  of  the  stochastic  controllability  matrix: 

Uk{ti)=  Qd(ti)  $(ti+i,ti)Qd(ti)  $^(t,-+i,ti)Q<i(^i)  ••• 

(2.55) 

As  the  eigenvalues  of  the  stochastic  controllability  matrix  become  large,  the  greater  the 
relative  benefit  of  smoothing. 

2.5  Digital  Filtering 

This  section  describes  the  types  of  digital  filters  considered  in  this  research:  Finite 
Impulse  Response  (FIR)  filters  and  Infinite  Impulse  Response  (HR)  filters. 

2.5.1  FIR  Filter.  The  causal  FIR  filter,  shown  in  Figure  2.2,  is  a  non-recursive 
discrete-time  filter  that  depends  only  on  the  present  and  past  input  signals.  The  transfer 
function  of  an  FIR  filter  is  given  as: 

JV-l 

Hiz)  =  Hn)  (2.56) 

n=0 
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where  h{n)  is  the  filter  impulse  response  of  length  N. 


Figure  2.2  Structure  of  FIR  Filter 

The  FIR  filter  has  several  interesting  properties: 

Linear  Phase  An  FIR  filter  can  be  designed  with  linear  phase,  wherein  we  assume  that 
there  will  be  no  phase  distortion  of  the  filtered  signal. 

Stability  FIR  filters  are  non-recursive  and  inherently  stable. 

Numerics  FIR  filters  are  less  sensitive  to  coefficient  accuracy  and  quantization  noise 
problems  than  IIR  filters. 

Filter  Order  To  achieve  the  same  magnitude  response  as  an  IIR  filter,  FIR  filters  are  of 
a  much  higher  order  than  IIR  filters. 

2.5.2  IIR  Filter.  The  IIR  filter,  shown  in  Figure  2.3,  is  a  recursive  discrete-time 
filter  that  uses  past  and  present  outputs  as  well  as  inputs.  The  transfer  function  of  an  IIR 
filter  is  given  as: 

H{z)  =  ' -  (2.57) 

The  use  of  past  and  present  outputs  in  the  IIR  filter  formulation,  also  known  as 
feedback,  results  in  some  interesting  properties: 

Non-Linear  Phase  A  casual  IIR  filter  has  a  non-linear  phase  response.  The  implemen¬ 
tations  of  a  non-causal  and  causal  IIR  filter  can  be  combined,  however,  to  produce  a 
zero  phase  response. 

Stability  IIR  filters  employ  feedback  and  can  be  unstable. 
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Figure  2.3  Structure  of  HR  Filter  (Canonic  Form) 


Numerics  IIR  filters  are  sensitive  to  coefficient  accuracy  and  quantization  noise  problems. 

Filter  Order  IIR  filters  are  typically  of  much  lower  order  than  a  corresponding  FIR  filter. 

Having  a  lower  order  also  gives  the  IIR  filter  a  better  (shorter)  transient  response, 

requiring  fewer  sample  periods  to  reach  steady-state  operation. 

There  are  four  classic  IIR  filters:  Butterworth,  Chebyshev,  inverse  Chebyshev  (or 
Chebyshev  II)  and  Elliptic.  These  four  filter  types  can  be  classified  under  two  error  ap¬ 
proximation  measures.  One  measure  involves  a  Taylor’s  series  expansion,  equating  as  many 
of  the  derivatives  of  the  desired  response  to  actual  response  as  possible.  The  other  measure 
attempts  to  minimize  the  maximum  difference  between  the  desired  and  actual  responses 
over  a  frequency  range.  These  IIR  filters  are  often  designed  as  analog  filters  and  trans¬ 
formed  to  digital  filters  using  the  bilinear  transformation.  Good  discussions  of  IIR  filter 
design  can  be  found  in  digital  signal  processing  textbooks,  such  as  [7,31]. 

2. 6  Summary 

This  chapter  has  given  a  brief  introduction  to  calculating  an  approximation  to  ve¬ 
locity  from  discrete  position  data.  Also  given  is  theory  relating  to  the  “raw”  GPS  mea¬ 
surements:  pseudorange,  deltarange  and  carrier- phase.  The  Kalman  filter /smoother  is 
presented  as  a  means  to  model  the  trajectory  of  an  aircraft  in  controlled  flight.  Finally, 
the  FIR  and  IIR  digital  filter  architectures  are  outlined  and  proposed  for  this  research.  The 
material  presented  in  this  chapter  forms  a  base  of  information,  essential  for  understanding 
the  results  to  be  presented  later. 


III.  Models 


3.1  Overview 

This  chapter  begins  by  presenting  the  simple  kinematic  models  describing  the  tra¬ 
jectory  of  an  aircraft  in  controlled  flight.  Next,  the  tropospheric  estimation  error,  receiver 
noise  and  multipath  error  models  are  developed  for  use  in  both  the  Kalman  filter/smoother 
and  discrete-time  filtering  and  numerical  methods  approaches.  Finally,  a  computer-assisted 
design  methodology  for  the  FIR  and  HR  discrete-time  filters  is  presented. 


3.2  Kinematic  Models 

The  kinematic  models  developed  in  this  section  are  intended  to  adequately  describe 
the  motion  of  an  aircraft  along  a  trajectory  segment.  A  trajectory  segment  is  defined  as  a 
continuous  portion  of  the  overall  trajectory. 

Kinematic  modeling  describes  motion  by  defining  a  kinematic  position  vector 

r(t)  =  {.■r(t),y(t),z(<)},  (3.1) 


a  kinematic  velocity  vector, 

Mi) 

^  dt  dt  '  dt  '  dt  ^ 


(3.2) 


a  kinematic  acceleration  vector, 

_  d?v{t)  _  d?x{t)  d?y{t)  d?z{t) 

’  dt^  dt^  ’  df^  ’  -f’ 


(3.3) 


and  so  on: 

d^rjt)  _  d»x(t)  dPy{t)  d^zjt) 
dt"  dt"  ’  dt"  ’  dt"  ^  ^ 

These  representations  of  motion  are  each  equivalent  if  the  appropriate  initial  conditions  are 

known  and  the  function  r(t)  is  continuously  differentiable.  In  practice,  these  assumptions 

generally  do  not  hold  and  the  closest  representation  to  the  measurement  is  chosen.  In  this 

way,  position  information  can  be  used  to  determine  an  approximation  to  velocity. 
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For  this  effort,  the  initial  conditions  are  not  precisely  known,  the  time  history  is 
discretized  and  the  position  measurements  are  noise  corrupted.  The  technique  of  pseudo¬ 
noise  addition  [25]  is  used  in  an  attempt  to  compensate  for  these  shortcomings.  The  models 
can  be  expressed  as  a  system  of  stochastic  difference  equations  of  the  form: 


x(i,)  =  -I-  Gd(t<)Wd(ti_i) 


(3.5) 


z{ti)  =  H(/i)x(ti)  +  v(ti). 


(3.6) 


where:  x(t,)  =  n-dimensional  system  state  vector 

$(/,-,  t,_i)  =  state  transition  matrix 

Gd(t,)  =  noise  distribution  matrix 

Wd(tt-i)  =  white  Gaussian  dynamics  noise  vector  of  covariance,  Q(t) 
z(t,)  =  m-dimensional  measurement  vector 

H(ti)  =  measurement  matrix 

v(ti)  =  white  Gaussian  measurement  noise  vector  of  covariance,  R(t). 
Three  simple  kinematic  models  are  used  in  this  research:  a  constant  velocity  with 
white  noise  added,  constant  acceleration  with  white  noise  added,  and  constant  acceleration 
with  noise  modeled  as  a  first  order  Markov  process. 


3.2.1  Constant  Velocity.  The  constant  velocity  model  assumes  that  the  aircraft’s 
acceleration  is  well  described  by  a  zero  mean  white  Gaussian  noise,  w^,  of  strength  Q^. 
The  constant  velocity  model  is  given  by: 


and 


X  = 


T 


Px  Py  Vy  Pz  V, 


^{ti,ti_i) 


1  T,  0  0  0  0 
0  1  0  0  0  0 
0  0  1  T,  0  0 
0  0  0  1  0  0 
0  0  0  0  1  r, 
0  0  0  0  0  1 


(3.7) 


(3.8) 
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where  is  the  sampling  time  of  the  discrete  system,  and  is  the  equivalent  discrete 
time  noise  approximated  by: 


and  Q(/)  is  given  by: 

qr  0  0 

Q(f)  =  0  0  (3.10) 

0  0 

where  q^,,  q^,  and  q^  are  tunable  parameters.  The  noise  distribution  matrix  is  given  by: 


G(t)  = 


0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

The  measurement  matrix  H  is  given  by: 


1  0  0  0  0  0 

H=  001000 
0  0  0  0  1  0 


(3.11) 


(3.12) 


with  a  measurement  noise  covariance  of: 


Tj;  0  0 

R  =  0  rj,  0 

0  0  r, 

where  Tj,,  Vy,  and  are  tunable  parameters. 


(3.13) 
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3.2.2  Constant  Acceleration.  Similarly  the  constant  acceleration  model  assumes 
jerk  is  modeled  as  a  zero  mean  white  Gaussian  noise  and  is  given  by: 

tT 

X  =  py  Vy  ay  p^  a,  (3-14) 

with  a  state  transition  matrix: 

IT^^OO  000  0 

OIT,  00  000  0 

00  100  000  0 

00  OIT,  ^00  0 

00  001  T,  00  0  (3.15) 

00  000  100  0 

00  000  01  r,  ^ 

00  000  001  r, 

00  000  000  1 

where  T,  is  the  sampling  time  of  the  discrete  system,  and  is  the  equivalent  discrete 
time  noise  given  by  (3.9)  and  (3.10).  The  noise  distribution  matrix  is: 

0  0  0 
0  0  0 
1  0  0 
0  0  0 

G(t)=  0  0  0.  (3.16) 

0  10 
0  0  0 
0  0  0 
0  0  1 


3-4 


The  measurement  matrix  H  is  given  by: 


100000000 

H=  000100000  (3.17) 

000000100 

with  a  measurement  noise  covariance  of: 

r*  0  0 

R=  0  rj,  0  (3.18) 

0  0  r, 

where  r^;,  r^,  and  are  tunable  parameters. 

S.S.S  Constant  Acceleration,  Noise  1st  Order  Gauss-Markov.  In  this  section, 
the  constant  acceleration  model  developed  in  Section  3.2.2  is  augmented  with  additional 
shaping  filter  states.  This  model  is  used  to  describe  the  jerk  of  the  aircraft  as  a  zero  mean, 
first  order  Gauss-Markov  process.  This  type  of  model  is  very  useful  for  approximating  a 
variety  of  empirically  observed  band-limited  noises  [24]  and  is  popular  in  the  open  literature 
[5,42,46].  The  model  is  described  by: 

x(t)  =  -(l/r)x(t)  +  w(t)  (3.19) 

Also  known  as  a  first-order  lag,  driven  by  white  Gaussian  noise  of  strength  Q  =  2rr^/r, 
this  process  produces  an  output  with  an  autocorrelation  described  by 

$(r)  =  E{x{t)x{t  A  r)}  =  (3.20) 

where  T  is  the  correlation  time  and  cr^  is  the  mean  squared  value.  The  first-order  Gauss- 
Markov  process  is  an  exponentially  time-correlated  process. 
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The  model  structure  is  identical  to  the  model  presented  in  Section  3.2.2  with  the 


exception  of  the  state  transition  matrix: 


where 


1  T, 
0  1 
0  0 


0  0 
0  0 
0  0 


0  0 
0  0 
0  0 


/l,3  =  /4,6  =  /7,9 
/2,3  =  /s.e  =  /s.g 
/3,3  =  /e.e  =  /9,9 


/i,3  0  0  0  0  0 

/2,3  0  0  0  0  0 

/3,3  0  0  0  0  0 

0  1  Z  /4,6  0  0 

0  0  1  /5,6  0  0 

0  0  0  /6,6  0  0 

0  0  0  0  1  T, 

0  0  0  0  0  1 

0  0  0  0  0  0 


T,  •  T  +  - 


0 

0 

0 

0 

0 

0 

/7,9 

/s,9 

/9,9 


1) 


(3.21) 


(3.22) 


3.3  Noise  Models 

The  major  GPS  error  sources  can  be  broken  down  into  three  categories:  Transmitter 
Errors,  Propagation  Errors,  and  Receiver  Errors.  Traditional  GPS  transmission  errors 
include  orbital  estimation  errors,  clock  errors  and  the  effects  of  Selective  Availability  (SA). 
Propagation  errors  include  the  effects  of  the  ionosphere  and  troposphere  on  the  GPS  signal, 
as  well  as  multipath  at  the  receiver  site.  Receiver  error  can  be  attributed  to  multipath 
effects  and  a  generic  receiver  noise  (Phase  Locked  Loop,  clock). 

For  the  pseudolite-based  positioning  system  used  in  this  research,  orbital  estimation 
and  SA  effects  are  not  applicable  sources  of  error.  Similarly,  since  the  pseudolite  signals 
do  not  travel  through  the  ionosphere,  this  source  of  error  is  neglected.  For  carrier-phase 
signals  the  pseudolite  clock  error  is  assumed  to  be  eliminated  by  the  differencing  approach 
in  Section  2.3.2. 1. 
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3.3.1  Tropospheric  Delay.  The  tropospheric  delay  effect  most  commonly  pre¬ 
sented  in  GPS  literature  actually  reflects  the  combined  effect  of  the  troposphere  and 
stratosphere.  The  troposphere  is  defined  as  the  first  10  km  above  the  Earth’s  surface. 
The  stratosphere  is  defined  as  the  region  just  above  the  troposphere  and  extends  to  60 
km  [19].  Figure  3.1  portrays  the  separation  of  the  troposphere  and  stratosphere. 


Since  most  of  the  water  vapor  is  concentrated  in  the  lower  atmosphere,  below  10 
km,  most  tropospheric  models  essentially  divide  the  tropospheric  effect  into  wet  and  dry 
effects  [16].  The  dry  component  accounts  for  80-90%  of  the  total  tropospheric  delay  effect 
and  is  predictable  within  1%  at  the  zenith  point,  while  the  wet  component  makes  up  10- 
20%  of  the  total  tropospheric  delay  effect  and  is  predictable  within  10-20%  at  the  zenith 
point  [19]. 

There  are  several  tropospheric  models  used  in  GPS  research  [16],  including  the  Saas- 
tamoinen,  Davis,  Goad  and  Goodman,  Yionoulis,  Black,  Chou,  Marini,  and  others.  Dif¬ 
ficulties  arise  in  using  the  traditional  GPS  tropospheric  models  in  this  research,  due  to 
the  inverted  nature  of  the  pseudolite  positioning  concept.  First,  the  GPS  tropospheric 
models  assume  that  the  GPS  signal  is  transmitted  from  satellites  in  high  orbits.  The  tro- 


3-7 


pospheric  delay  is  then  calculated  from  the  top  of  the  stratosphere  to  the  receiver;  in  this 
case,  an  aircraft  flying  somewhere  in  the  troposphere.  However,  the  delay  term  that  we 
are  concerned  with  is  between  the  aircraft  and  a  ground  receiver.  One  might  compensate 
for  this  by  calculating  the  total  delay  between  the  satellite  and  ground  receiver  and  then 
subtracting  the  delay  between  the  satellite  and  aircraft.  Secondly,  the  tropospheric  delay 
is  calculated  as  a  function  of  the  zenith  angle.  As  these  zenith  angles  go  beyond  70  —  80° 
many  of  the  GPS  tropospheric  models  attempt  to  model  the  curvature  of  the  Earth  and 
consequently  overestimate  the  delay  term  [16]. 

3.3.2  Receiver  Noise  and  Multipath  Error.  An  approximation  to  the  strength  of 
receiver  noise  is  to  treat  1%  of  the  signal  wavelength  as  white  noise  [19].  This  noise  is  due 
to  the  processing  of  the  signals  by  the  receiver  hardware.  For  example,  the  GPS  LI  carrier 
is  transmitted  at  1573.42  MHz,  corresponding  to  a  wavelength  of  19cm.  Thus,  the  receiver 
noise  in  the  carrier- phase  Ll  measurement  can  be  modeled  as  white  noise  with  a  strength 
of  1.9mm.  Similarly,  the  code  phase  frequency  is  1.023  MHz  resulting  in  a  wavelength  of 
293m  and  a  receiver  noise  strength  of  2.93m. 

Multipath  effect  is  caused  by  the  antenna  being  exposed  to  reflected  signals,  causing 
interference  in  the  receiver  [19].  The  effect  of  multipath  is  highly  dependent  upon  the 
geographic  and  cultural  features  of  the  receiver  site.  Several  parameters  are  used  to  describe 
the  reflection  of  a  signal  off  the  Earth’s  surface.  These  include  the  dielectric  constant 
(A’r^),  conductivity  (ct)  and,  index  of  refraction  (n  =  y/Krg).  The  multipath  effect  is 
also  dependent  upon  the  elevation  angle  of  the  pseudolite,  0^.  Figure  3.2  shows  the  effect 
of  reflection.  The  critical  elevation  angle  is  defined  as  the  angle  below  which  Snell’s  law 
predicts  reflection.  Using  Snell’s  law,  the  critical  angle  [9^)  is  determined  as: 

Oa  =  sin"^(— )  (3.23) 

Ug 

6,  =  90°  -  9^  (3.24) 

where  n-a  is  the  index  of  refraction  of  air  and  Hg  is  the  index  of  refraction  of  the  ground. 
Sandy,  dry  areas,  typical  of  the  desert  of  New  Mexico,  have  an  index  of  refraction  Ug  =  3.2 
[19].  The  index  of  refraction  of  air  is  defined  as  Ua  =  1.  Using  (3.24)  and  (3.24),  the  critical 
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Figure  3.2  Multipath  Signal  Reflection 


elevation  angle  is  determined  to  be  6^  =  71. S”.  Thus,  as  aircraft  flying  over  the  RRS  range 
become  closer  to  the  ground  and  fly  farther  away  from  the  receivers,  they  will  induce 
multipath  effects  at  elevation  angles  less  than  Oc-  Section  4.2  illustrates  the  relationship 
between  the  flight  profile,  receiver  locations,  and  geometry  concerns  (e.g.  elevation  angles). 

3.3.3  Noise  Model  Implementation.  Figure  3.3  is  a  SIMULINK  block  diagram 
of  the  noise  model  used  in  this  research.  This  model  outputs  the  ranging  error  between 
a  given  pseudolite  transmitter  and  receiver.  Inputs  into  the  model  are  a  time  history  of 
range  and  elevation  angle. 

There  are  three  sections  of  the  noise  model:  the  multipath  error,  receiver  noise  and 
tropospheric  estimation  error.  The  multipath  error  and  tropospheric  estimation  error  are 
modeled  in  this  research  as  first  order  Markov  processes  with  time  constants  of  5  seconds 
and  30  seconds  respectively.  Receiver  noise  is  customarily  modeled  as  “white”  noise.  The 
multipath  error  model  uses  the  elevation  angle  as  a  parameter  to  scales  the  multipath 
range  error  with  the  following  formula: 

I  1  —  sin(|0ei))  when  Be  <  71.8° 

1  0,  when  9^  >  90° 
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Figure  3.3  Block  Diagram  of  Noise  Model 
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Similarly,  the  tropospheric  estimation  error  model  takes  range  as  a  parameter  and  scales 

^irop  • 

■  lOfcm  ^  ^ 

such  that  the  error  is  normalized  about  the  radius  of  the  troposphere,  or  10  km. 


Two  noise  models  are  used  in  this  research.  The  first  model,  hereafter  refered  to  as 
the  “colored”  noise  model,  contains  tropospheric,  multipath  and  receiver  noise  error  and 
is  described  above.  The  other  model  is  a  “white”  noise  model,  where  the  tropospheric 
and  multipath  effects  are  neglected.  The  receiver  noise  strength  is  increased  to  achieve  an 
overall  0.1m  3-D  RMS  position  solution. 

Fifty  realizations  of  each  noise  model  were  generated  and  saved.  This  allows  both 
the  digital  filtering  and  numerical  differentiation  approach  and  Kalman  filtering/smoothing 
approache  to  share  the  same  noise  realizations. 


3.4  Digital  Filter  Design 

In  this  research,  the  digital  filters  were  synthesized  with  the  use  of  computer-aided 
design  tools  found  in  MATLAB’s  Signal  Processing  Toolbox.  The  design  procedure  involves 
identifying  the  following  parameters  that  characterize  the  performance  of  a  low  pass  filter: 

•  Rp  =  allowed  passband  loss 

•  i?j,  =  desired  stopband  attenuation 

•  6p  =  allowed  passband  ripple 

•  Sg  =  allowed  stopband  ripple 

•  fp  =  desired  passband  frequency 

•  /s  =  desired  stopband  frequency 

These  design  parameters  are  shown  in  Figure  3.4  The  passband  includes  all  positive  fre¬ 
quencies  up  to  fp ,  the  stopband  includes  all  frequencies  above  fg ,  and  the  transition  band 
is  defined  as  the  region  between  fp  and  /,.  For  FIR  filters,  the  MATLAB  function  remezord 
accepts  as  inputs:  fp,  fg,  6p  and  6,  and  returns  the  parameters  required  to  design  the  FIR 
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Figure  3.4  Parameters  for  Filter  Design 


filter  using  the  well-known  Parks-McClellan  algorithm.  Similarly,  HR  filters  are  designed 
using  MATLAB  functions  that  are  capable  of  generating  Butterworth,  Chebyshev  (I  and 
II)  and  Elliptic  filters.  By  varying  combinations  of  the  parameters  Rp,  R,,  fp,  fg,  6p  and 
6s,  the  filter  design  is  accomplished  in  an  iterative  and  intuitive  manner. 

For  example,  to  design  the  lowest  order  Chebyshev-I  low  pass  filter  whose  passband 
frequency  is  0.5,  with  stopband  frequency  of  0.8,  with  no  more  than  1  dB  of  loss  in  the 
passband,  has  at  least  20  dB  attenuation  in  the  stopband,  and  allows  0.1  dB  of  passband 
ripple,  one  would  use  the  MATLAB  commands  cheblord  and  chebyl,  as  in  the  following: 

[M,  Wn]  =  cheblord(0.5,  0.8,  1,  20); 

[B,  A]  =  chebyl (N,  Wn,  0.1); 

Now  the  filter  is  completely  described  by  vectors  B  and  A,  representing  the  numerator  and 
denominator  coefficients  of  the  filter  transfer  function. 
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IV.  Results 


4.1  Overview 

This  chapter  presents  the  results  of  the  development  and  testing  of  the  two  ap¬ 
proaches  taken  in  this  research:  Digital  filtering  with  numeric  differentiation  and  Kalman 
filtering/smoothing.  The  chapter  begins  with  a  discussion  of  the  aircraft  trajectories  used 
in  this  research. 

4.2  Aircraft  Trajectory 

The  software  package  PROFGEN  [32]  is  used  to  create  the  flight  profile  used  in 
this  research.  The  flight  profile  is  essentially  the  same  used  in  previous  CIGTF-sponsored 
AFIT  thesis  research.  It  is  intended  to  simulate  a  two-hour  mission  performed  by  a  high 
performance  multi-role  fighter  aircraft.  Figure  4.1  shows  a  two  dimensional  representation 
of  the  profile.  Figure  4.2  shows  the  time  history  of  “gees”  as  experienced  by  the  aircraft, 
which  provides  an  overall  measure  of  the  dynamics  of  the  profile. 


Figure  4.1  2-D  View  of  Flight  Profile 


Figure  4.3  shows  a  three-dimensional  view  of  the  profile,  which  is  oriented  over  the 
CIGTF  RRS  range.  Overlayed  on  the  profile  are  the  locations  of  the  six  receiver  sites 
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considered  for  this  research.  Table  4.1  lists  the  geographic  location  of  the  receivers.  These 
receivers  are  located  at  existing  RRS  transponder  sites  [33],  and  are  not  formally  optimized 
for  geometry.  The  position  of  the  receivers  relative  to  the  flight  profile  can  also  be  seen  in 
Figure  4.4. 


Table  4.1  Pseudolite  Receiver  Locations  for  NRS  Simulation 


# 

Location 

Longitude 

Latitude 

Altitude 

R1 

Tula  PK,  NM 

33'’01.36' 

-106‘’08.20' 

1322.53m 

R2 

TDC,  NM 

32‘’55.58' 

-106“08.50' 

1241.76m 

R3 

Oscura  PK,  NM 

33M4.58' 

-106‘’22.14' 

2417.51m 

R4 

Salinas,  NM 

33n7.55' 

-106“31.44' 

2695.11m 

R5 

Sac  Peak,  NM 

32'’47.16' 

-105M9.15' 

2804.81m 

R6 

Twin  Buttes,  NM 

32‘’42.12' 

-106‘’07.38' 

1365.71m 

Figure  4.4  Top  View  of  Flight  Profile  and  Receiver  Locations 

Of  concern  is  the  need  to  keep  the  configuration  of  the  receivers  non-coplanar,  to 
ensure  a  reasonable  three  dimensional  solution  [36].  Furthermore,  the  ranges  and  elevation 
angles  between  the  receivers  and  the  aircraft  are  of  interest,  as  they  effect  the  tropospheric 
delay  and  multipath  errors.  Figure  4.5  shows  the  ranges  and  elevation  angles  for  the  flight 
profile.  Notice  that  the  elevation  angles  are  typically  below  10  degrees,  with  the  exception 
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of  when  the  aircraft  flies  almost  directly  over  the  receivers.  These  low  elevation  angles  are 
not  only  good  conditions  for  inducing  multipath  errors,  but  also  represent  an  undesirable 
coplanar  receiver  configuration. 

There  is  a  segment  of  the  trajectory  that  shows  small  negative  elevation  angles. 
Elevation  angle  is  calculated  using  a  plane  tangent  to  the  surface  of  the  Earth  as  reference. 
Since  altitude  is  computed  as  height  above  the  WGS-84  reference  ellipsoid,  these  small 
negative  angles  are  a  result  of  the  aircraft  altitude  being  below  the  tangent  plane  extending 
out  over  a  curved  Earth. 


Figure  4.5  Range  and  Elevation  Angles  Between  Aircraft  and  Receivers 

Another  metric  for  the  geometry  of  receiver  locations  is  to  consider  the  Geometric 
Dilution  of  Precision  (GDOP)  .  GDOP  is  a  ratio  of  3-D  position  (and  time)  accuracy  to 
measurement  accuracy  [18]. 

GDOP  is  defined  as: 

GDOP  =  yTr{(AT  .  A)-i}  (4.1) 
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where  the  matrix  A  is  defined  as, 


Dxo 

Dyo 

Dzq 

Dt^ 

Dxi 

Dyi 

Dzi 

Dt 

-Dxyj  Dyfi  DZfi  Dtfi 


(4.2) 


containing  the  directional  derivatives: 


Dxi  = 

{Xj  -  Xi) 

(4.3) 

-  Xif  -t-  {yj  -  Vif  +  {Zj  -  Zif 

Dyi  = 

(y,  -  yi) 

(4.4) 

ys” 

-  +  (yj  -  y>Y  +  -  ^iY 

Dzi  - 

(Zj  -  Zi) 

(4.5) 

-  3:iY  +  (%•  -  yiY  +  -  ^iY 

Dti  = 

1 

(4.6) 

GDOP  is  calculated  for  the  entire  flight  profile  using  the  receiver  locations  in  Table 
4.1  and  is  presented  in  Figure  4.6.  As  GDOP  groves,  the  set  of  equations  that  must  be 
solved  simultaneously  to  yield  position  and  time  become  more  linearly  dependent.  Not  only 
does  the  uncertainty  in  the  solution  grow  as  GDOP  becomes  large,  but  in  the  presence 
of  measurement  noise,  the  existence  of  a  solution  is  in  jeopardy.  In  practice,  acceptable 
GDOP  for  a  GPS  receiver  is  considered  to  be  an  arbitrary  value  usually  below  10.  Figure 
4.7  shows  a  portion  of  the  flight  segment  with  reasonable  GDOP. 

In  simulation  it  is  possible  to  preserve  a  good  GDOP  throughout  the  flight  by  placing 
additional  receivers  along  the  flight  path.  A  crude,  yet  effective  algorithm  was  developed 
to  accomplish  this.  The  algorithm  monitors  GDOP  during  a  flight.  If  GDOP  exceeds  a 
given  threshold,  a  search  in  a  fixed  radius  around  the  aircraft  locates  the  site  yielding  the 
best  improvement  in  GDOP.  The  number  of  receivers  used  to  determine  GDOP  at  one 
time  is  fixed,  and  the  addition  of  a  new  receiver  generally  results  in  the  removal  of  an 
existing  receiver.  As  receiver  sites  exceed  range  or  elevation  limits,  they  are  removed  from 
consideration  as  well. 
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Figure  4.7  Plan  View  of  Flight  Profile  and  GDOP  for  Truncated  Flight  Segment 
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Figure  4.8  Plan  View  of  Receiver  Locations  Ensuring  GDOP  <  10 

In  Figure  4.8  we  see  the  receiver  laydown  pattern  for  the  entire  flight  profile.  During 
this  flight,  receivers  were  placed  such  that  GDOP  was  not  allowed  to  exceed  10,  as  shown 
in  Figure  4.9.  The  receiver  locations  are  marked  with  an  X.  The  pattern  in  Figure  4.8, 
forms  a  “racetrack”  surrounding  the  flight  profile.  The  relative  proximity  of  the  receiver 
locations  to  the  flight  profile  suggests  that  rather  than  trying  to  encompass  the  entire  flight 
profile  with  a  single  set  of  receivers,  a  “racetrack”  pattern  may  be  a  potential  strategy  for 
SARS  receiver  deployment^.  This  algorithm  is  effective  in  limiting  GDOP  for  simulation 
purposes  yet  also  provides  insight  into  the  real-world  geometry  problem  facing  SARS. 

4-3  Digital  Filtering  with  Numeric  Differentiation 

In  this  section,  the  validation  and  results  of  the  digital  filtering  with  numeric  differ¬ 
entiation  approach  are  presented.  First  each  of  the  components  of  the  block  diagram  are 
validated  for  performance.  The  validation  procedure  involves  benchmarking  the  compo- 

^The  issue  of  receiver  geometry  for  a  realistic  flight  profile  and  a  reasonable  number  of  receiver  sites 
should  be  addressed  in  more  detail  as  an  issue  of  feasibility  for  the  SARS  project. 
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Figure  4.9  GDOP  for  Entire  Flight  Profile  Using  Receiver  Locations  in  Figure  4.8 

nents  in  the  absence  of  noise,  in  the  presence  of  white  noise,  and  finally  in  the  presence  of 
a  more  realistic,  colored  noise. 

4.3.1  Analysis  of  Sampling  Rates  for  Alias  Effects.  Since  the  time  history  of 
GPS  position  data  is  treated  as  a  sampled  signal,  it  is  necessary  to  ensure  that  the  Nyquist 
sampling  criterion  is  met.  A  short  flight  profile  was  generated  that  simulated  the  highest 
anticipated  flight  dynamics;  minimum  radius  9  g  s-turns.  This  flight  profile  was  generated 
in  ECEF  coordinates  at  two  sampling  rates:  1000  Hz  and  10  Hz.  The  discrete  Fourier 
transform  was  implemented  via  MATLAB’s  built-in  Fast  Fourier  Transform  (FFT)  routine 
to  examine  the  frequency  spectra  of  the  signal.  The  flight  profile  was  generated  at  1000 
Hz  to  ensure  that  if  any  aliasing  was  present  at  lower  sampling  rates,  it  could  be  detected. 

Figures  4.10  and  4.11  show  the  frequency  spectra  for  1000  Hz  and  10  Hz  respectively. 
Figure  4.10  verifies  that  there  are  no  components  of  the  frequency  spectra  at  high  frequen¬ 
cies.  Similarly  at  10  Hz,  there  is  no  evidence  of  frequency  components  that  might  suggest 
aliasing  at  lower  frequencies. 
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4-3.2  Numerical  Differentiation  Algorithm  Performance.  Two  procedures  for 
determining  a  numerical  derivative  were  examined.  The  first  procedure  involves  an  ex¬ 
trapolated  Taylor’s  series  expansion  of  the  derivative  operation  as  described  in  Section 
2.2.2.  The  second  procedure  uses  splines  to  perform  a  curve  fit  to  the  numerical  data.  An 
analytical  derivative  of  the  equation  of  the  splines,  evaluated  at  the  points  of  interest,  gives 
an  approximation  to  the  derivative. 

4.3.2. 1  Central  Difference  Equations.  The  central  difference  equation  ap¬ 
proximation  for  the  first  derivative  is  well  documented  [6,15,22,45]  in  numerical  methods 
texts.  Equations  (2.10),  (2.11),  and  (2.12)  show  the  algorithms  used  in  this  research. 

To  test  the  algorithms,  a  test  signal  is  defined: 

r(x)  =  sin(a;^),  (4.7) 

which  has  an  exact  analytic  derivative: 

r'(a;)  =  2a;  cos(a:^).  (4-8) 

Equations  (4.7)  and  (4.8)  are  evaluated  over  the  interval  [0,37r]  at  a  given  sample  rate. 
The  numerical  derivative  is  computed  and  is  compared  to  the  exact  value  determined  from 
(4.8).  As  expected,  (2.12)  demonstrates  better  performance  than  the  other  equations  in 
the  absence  of  noise.  Figure  4.12  shows  the  error  committed  by  the  first,  second  and  third 
order  central  difference  equations  at  a  100  Hz  sample  rate.  We  can  see  from  Figure  4.12 
that  the  third  order  central  difference  equation  shows  small  errors  on  the  order  of  10“®, 
compared  to  10“^  and  10“^  for  the  first  and  second  order  central  difference  equations 
respectively. 

The  test  signal  is  chosen  such  that  as  x  increases,  the  frequency  of  the  test  signal  si¬ 
nusoid  increases,  and  the  approximation  to  the  derivative  worsens.  This  error  is  attributed 
to  two  factors:  the  truncation  of  the  higher  order  terms  in  the  Taylor’s  series  and  the 
nature  of  the  sampled  waveform.  Truncation  establishes  the  order  of  magnitude  of  the 
error,  while  the  increasing  frequency  of  the  signal  exposes  the  error  due  to  sampling. 
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Figure  4.12  Error  in  Central  Difference  Approximations  in  the  Absence  of  Noise 
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Next  we  examine  the  performance  of  these  algorithms  under  realistic  dynamics,  as 
provided  by  the  PROFGEN  flight  profile.  The  PROFGEN  indicated  positions  are  numer¬ 
ically  differentiated  using  the  central  difference  equations.  These  computed  velocities  are 
compared  to  the  PROFGEN  indicated  velocities,  which  are  considered  to  be  truth  data. 
See  Appendix  A  for  a  description  of  the  error  metrics  used  in  this  research.  These  results 
are  collected  in  Table  4.2.  As  expected,  the  higher  order  algorithm  provides  the  best  per¬ 
formance  at  both  data  rates.  Also,  as  the  data  rate  is  increased  from  1  Hz  to  10  Hz,  all  of 
the  numerical  differentiation  algorithms  improve. 


Direction 

Rate 

Error  (m/s) 

Icr  (m/s) 

Order 

3-D  RMS 

1  Hz 

6.125160e-02 

4.035006e-01 

3-D  RMS 

1  Hz 

1.355792e-02 

1.164690e-01 

3-D  RMS 

1  Hz 

8.047649e-03 

8.551035e-02 

3-D  RMS 

6.510429e-04 

4.601638e-03 

3-D  RMS 

ill 

5.020955e-06 

1.171970e-04 

3-D  RMS 

Hi 

1.301999e-06 

5.134946e-05 

Table  4.2  Performance  of  Numeric  Differentiators  Using  Flight  Profile 


Figure  4.13  shows  the  3-D  RMS  error  of  the  third  order  numerical  differentiator  for  a 
10  Hz  data  rate.  This  plot  shows  large  “spikes”  of  error  at  specific  time  intervals,  which  are 
correlated  in  time  to  the  highly  dynamic  portions  of  the  flight  trajectory  as  seen  in  Figure 
4.2.  While  the  magnitudes  of  the  spikes  in  Figure  4.13  vary  for  the  different  numerical 
differentiation  algorithms  and  data  rates,  the  location  of  these  “spikes”  remains  consistent. 

We  would  then  expect  better  performance  of  the  numerical  differentiation  routines 
during  periods  of  low  dynamics.  Table  4.3  shows  the  numeric  differentiator  performance 
over  a  segment  of  the  flight  profile  with  benign  dynamics  (t  =  0  to  100  seconds).  Comparing 
Tables  4.2  and  4.3  we  see  a  93%-98%  improvement  in  3-D  RMS  error. 

While  it  is  worthwhile  to  consider  the  performance  of  these  numerical  differentiators 
in  the  absence  of  noise,  we  can  also  determine  the  performance  of  each  algorithm  in  the 
presence  of  corrupting  noise.  Two  types  of  noise  corruption  are  employed.  One  is  a 
simulation  of  “white”  Gaussian  noise,  while  the  other  is  a  combination  of  white  and  time 
correlated  noise,  hereafter  referred  to  as  “colored”  Gaussian  noise.  The  each  noise  type  is 
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Figure  4.13  Error  in  3rd  Order  Numerical  Derivative  Using  Flight  Profile  at  10  Hz 


Direction 

Rate 

Mean  Error  (m/s) 

1(7  (m/s) 

Order 

3-D  RMS 

1  Hz 

4.106344e-03 

1.386701e-02 

3-D  RMS 

1  Hz 

2.907217e-04 

3.474643e-03 

3-D  RMS 

1  Hz 

1.330864e-04 

2.321516e-03 

3-D  RMS 

10  Hz 

4.145871e-05 

1.418019e-04 

1st 

3-D  RMS 

10  Hz 

2.859088e-07 

1.245242e-05 

2nd 

3-D  RMS 

10  Hz 

1.123711e-07 

7.592645e-06 

3rd 

Table  4.3  Performance  of  Numeric  Differentiators  Using  Benign  Flight  Profile 
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tuned  to  provide  an  overall  0.1  m  3-D  RMS  lo  error  in  position  as  described  in  Section 
3.3.3. 


Direction 

Rate 

Mean  Error 
Colored  (m/s) 

la  (m/s) 

Mean  Error 
White  (m/s) 

la  (m/s) 

Order 

3-D  RMS 

1  Hz 

4.033391e-01 

9.760447e-03 

4.568517e-01 

3.114718e-02 

3-D  RMS 

1  Hz 

5.264175e-01 

1.240445e-02 

5.284405e-01 

4.333967e-02 

3-D  RMS 

1  Hz 

5.949750e-01 

1.387963e-02 

5.467742e-01 

5.120870e-02 

3-D  RMS 

10  Hz 

6.300444e-02 

4.476222e-03 

9.252125e-02 

8.219554e-01 

1st 

3-D  RMS 

10  Hz 

6.256138e-02 

5.565853e-03 

1.231328e-01 

1.114507e-f00 

2nd 

3-D  RMS 

10  Hz 

6.930055e-02 

5.707011e-03 

1.449810e-01 

1.306353e-}-00 

3rd 

Table  4.4  Performance  of  Numeric  Differentiators  in  the  Presence  of  Noise  Using  Benign 
Flight  Profile 


The  results  in  Table  4.4  show  the  performance  of  the  three  central  difference  equation 
algorithms  in  the  presence  of  noise,  for  a  50  Monte  Carlo  run  simulation.  For  both  the  1 
Hz  and  10  Hz  data  rates  the  benefit  of  using  the  higher  order  differentiators  is  no  longer 
evident.  In  fact,  at  1  Hz  the  higher  order  differentiators  show  a  14%  to  47%  increase  in  3-D 
RMS  error  compared  to  the  first  order  equation.  At  10  Hz,  the  second  order  differentiator 
slightly  outperforms  the  first  order  algorithm  for  the  case  of  colored  noise.  However,  in  all 
other  cases  at  10  Hz,  the  first  order  differentiator  provides  the  best  performance. 

4. 3. 2. 2  Differentiation  Using  Cubic  Splines.  Another  method  of  approxi¬ 
mating  a  derivative  is  to  determine  a  closed  form  expression  for  a  curve  that  fits  the  data 
to  be  operated  on.  An  analytic  derivative  is  taken  of  the  curve  and  then  evaluated  at  the 
points  of  interest  to  generate  the  derivative  approximation.  While  there  are  many  curve 
fitting  methods  available,  cubic  splines  handle  large  sets  of  data  points  well  and  guarantees 
continuous  derivatives  over  the  entire  data  set  [15,45].  As  in  the  previous  section,  (4.7)  and 
(4.8)  are  evaluated  for  a:  =  0  to  37r  at  a  given  sample  rate.  The  derivative  is  approximated 
using  cubic  splines  and  is  compared  to  the  exact  value  determined  from  (4.8).  Comparing 
Figures  4.12  and  4.14  shows  the  cubic  spline  derivative  performing  between  the  second  and 
third  order  numerical  differentiators. 
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Sample  Number 

Figure  4.14  Error  in  Cubic  Spline  Derivative  in  the  Absence  of  Noise 

As  with  the  numerical  dilferentiators,  it  is  of  interest  to  examine  the  performance 
of  the  cubic  spline  method  under  a  realistic  flight  profile.  The  cubic  spline  method  shows 
significant  improvement  in  velocity  accuracy  as  the  data  rate  increases  from  1  Hz  to  10 
Hz.  Also,  under  this  realistic  flight  profile,  the  cubic  spline  method  slightly  outperforms 
the  third  order  numerical  differentiator.  The  results  of  the  cubic  spline  differentiator  are 
tabulated  in  Table  4.5.  In  Figure  4.15  the  3-D  RMS  error  of  the  cubic  spline  differentiator 
for  10  Hz  is  shown.  The  location  and  magnitude  of  the  large  spikes  compares  favorably  to 
what  is  shown  in  Figure  4.13  as  well  as  the  aircraft  dynamics  shown  in  Figure  4.2. 

We  can  also  compare  the  performance  of  the  cubic  spline  differentiator  to  the  numer¬ 
ical  differentiation  routine,  within  the  benign  portion  of  the  flight  profile.  The  performance 
advantage  of  the  cubic  spline  method  is  easily  seen  by  comparing  Tables  4.3  and  4.5.  The 
cubic  spline  differentiator  shows  an  improvement  over  the  third  order  central  difference 
equation  in  3-D  RMS  error  of  64%  and  84%  for  1  Hz  and  10  Hz  data  rates  respectfully. 

Cubic  spline  differentiation  can  also  be  applied  to  noisy  signals.  To  accommodate 
noise,  a  smoothing  parameter,  p  £  [0---1],  is  introduced.  For  p  =  0,  the  cubic  spline 
algorithm  performs  a  least-squares  straight  line  fit  to  the  data.  For  p  =  1,  the  natural 
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Figure  4.15  Error  in  Cubic  Spline  Derivative  Using  Flight  Profile  at  10  Hz 


Direction 

Rate 

Mean  Error  (m/s) 

1(7  (m/s) 

Profile 

3-D  RMS 
3-D  RMS 

1  Hz 

10  Hz 

7.373304e-03 

1.345184e-06 

7.608367e-02 

4.488712e-05 

3-D  RMS 
3-D  RMS 

1  Hz 

10  Hz 

4.702409e-05 

1.767529e-08 

4.779469e-04 

1.428987e-07 

Truncated 

Truncated 

Table  4.5  Performance  of  Cubic  Spline  Differentiator  in  the  Absence  of  Noise 
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cubic  spline  fit  is  performed.  As  p  moves  from  0  to  1,  the  cubic  spline  algorithm  changes 
from  one  extreme  to  the  other.  The  smoothing  parameter  is  in  general  arbitrarily  chosen. 
In  this  case,  p  =  0.85  yielded  the  best  performance. 


Direction 

Rate 

Mean  Error 
Colored  (m/s) 

1<T  (m/s) 

Mean  Error 
White  (m/s) 

Icr  (m/s) 

Smoothing 

Parameter 

3-D  RMS 

1  Hz 

5.241575e-02 

4.872613e-03 

9.304974e-02 

7.544461e-03 

3-D  RMS 

10  Hz 

8.828444e-02 

6.173212e-03 

8.130258e-02 

5.698531e-03 

Table  4.6  Performance  of  Cubic  Smoothing  Spline  Differentiators  in  the  Presence  of  Noise 
Using  Benign  Flight  Profile 


Table  4.6  shows  the  performance  of  the  cubic  spline  differentiator  in  the  presence  of 
noise.  Comparing  the  cubic  spline  differentiator  results  to  the  central  difference  equation 
differentiator  results  of  Table  4.4  favors  the  cubic  spline  differentiator.  Except  for  the  case 
of  colored  noise  at  10  Hz,  the  cubic  spline  yields  lower  3-D  RMS  errors  by  40%  to  80%. 

4-3.3  Digital  Filter  Results.  In  order  to  validate  the  performance  of  the  digital 
filters  used  in  this  research,  we  first  examine  what  penalty,  if  any,  results  from  their  use. 
For  example,  filter  cutoff  frequencies  must  be  chosen  to  eliminate  as  much  noise  as  possible, 
while  preserving  the  true  underlying  signal.  To  do  this,  we  apply  low  pass  filters  to  both 
noisy  and  noiseless  signals,  and  observe  the  results.  In  both  cases,  the  difference  between 
the  original,  noiseless  signal  and  the  filtered  signal  is  the  error  we  wish  to  quantify. 

4-3.3. 1  IIR  Filter  Performance.  Four  types  of  HR  filters  were  examined: 
Chebyshev  I,  Chebyshev  II,  Butterworth  and  Elliptic.  The  filters  were  evaluated  for  per¬ 
formance  at  data  rates  of  1  Hz  and  10  Hz  using  the  truncated  flight  profile.  The  filter 
design  parameters  were  adjusted  by  trial  and  error  for  each  filter  type  to  achieve  it’s 
“best”  performance. 

Table  4.7  lists  the  design  parameters  and  resulting  filter  order  used  for  each  of  the 
filter  designs.  The  passband  frequencies  (/p)  and  stopband  frequencies  (fs)  are  given  in 
normalized  frequency  (0  to  1),  the  allowable  passband  ripple  {6p)  and  desired  stopband 
attenuation  (6*)  are  given  in  decibels.  These  design  parameters  are  used  to  describe  the 
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desired  properties  of  a  low  pass  filter,  which  is  synthesized  using  MATLAB  as  described 
in  Section  3.4. 


Type 

RpidB) 

RsidB) 

fp 

fs 

6p  (dB) 

<5.  (dB) 

order 

Chebyshev-I 

QRgjjjl 

0.55 

IBi 

le-07 

- 

9 

Chebyshev-II 

0.55 

- 

le-07 

8 

Butterworth^ 

le-07 

30 

0.55 

0.8 

- 

- 

13 

Elliptical 

le-07 

40 

0.55 

0.8 

le-07 

- 

7 

Table  4.7  HR  Filter  Design  Parameters  and  Resulting  Filter  Order 


The  Butterworth  filter  designed  using  parameters  similar  to  the  other  HR  filters 
was  unstable.  Since  the  Chebyshev  (I  and  II)  and  Elliptic  filters  would  provide  better 
performance  than  a  reduced  order  Butterworth,  the  Butterworth  filter  was  removed  from 
consideration.  Table  4.8  shows  the  results  of  applying  the  three  types  of  filters  to  noiseless 
position  data  at  1  Hz  and  10  Hz  rates  respectively.  The  Elliptic  filter  trails  the  performance 
of  the  Chebyshev  filters,  with  the  Chebyshev-II  filter  yielding  the  best  performance. 


Direction 

Filter  Type 

Rate 

Mean  Error  (m) 

3-D  RMS 

Chebyshev-I 

4.112376e-02 

3-D  RMS 

Chebyshev-II 

1.058158e-03 

3-D  RMS 

Elliptic 

9.470536e-02 

3-D  RMS 

Chebyshev-I 

10  Hz 

3.677020e-04 

3-D  RMS 

Chebyshev-II 

10  Hz 

1.521616e-04 

3-D  RMS 

Elliptic 

10  Hz 

8.528341e-04 

Table  4.8  Performance  of  HR  Filters  in  the  Absence  of  Noise 


The  performance  of  the  HR  filters  in  the  presence  of  noise  is  displayed  in  Table  4.9. 
Tests  were  performed  at  1  Hz  and  10  Hz  data  rates  with  both  white  noise  and  colored  noise 
models.  Although  the  Chebyshev-II  filter  gave  the  best  results  in  the  absence  of  noise,  the 
Chebyshev-I  filter  outperformed  all  other  filters  except  for  the  case  of  colored  noise  at  1  Hz. 
The  fact  remains,  however,  that  the  prefilter  is  unable  to  significantly  reduce  the  amount 
of  noise  in  the  position  signal  before  numerical  differentiation  takes  place. 


^Unstable  filter 


Direction 

Filter  Type 

Rate 

Mean  Error  (m) 

Noise 

3-D  RMS 

Chebyshev-I 

na 

1.387863e-01 

White 

3-D  RMS 

Chebyshev-II 

1.432184e-01 

White 

3-D  RMS 

Elliptic 

na 

1.810501e-01 

White 

3-D  RMS 

Chebyshev-I 

10  Hz 

1.173286e-01 

3-D  RMS 

Chebyshev-II 

10  Hz 

1.368191e-01 

3-D  RMS 

Elliptic 

10  Hz 

1.190525e-01 

White 

3-D  RMS 

Chebyshev-I 

■EB 

1.883719e-01 

Colored 

3-D  RMS 

Chebyshev-II 

1.692749e-01 

Colored 

3-D  RMS 

Elliptic 

lli 

2.291981e-01 

Colored 

3-D  RMS 

Chebyshev-I 

10  Hz 

1.384984e-01 

Colored 

3-D  RMS 

Chebyshev-II 

10  Hz 

1.413331e-01 

Colored 

3-D  RMS 

Elliptic 

10  Hz 

1.390585e-01 

Colored 

Table  4.9  Performance  of  HR  Filters  in  the  Presence  of  Noise 


4. 3. 3. 2  FIR  Filter  Performance.  While  the  HR  prefilter  was  able  to  reduce 
the  overall  error,  albeit  slightly,  an  FIR  filter  could  not  be  designed  to  meet  the  performance 
requirements.  The  MATLAB  routine  remezord  was  used  to  implement  an  FIR  filter  with 
the  characteristics  seen  in  Table  4.10.  Comparing  Tables  4.7  and  4.10,  the  order  of  the  FIR 


Type 

Rp{dB) 

Rs{dB) 

fp 

fs 

EnsEa 

6s  (dB) 

order 

FIR 

- 

IQI 

- 

68 

FIR 

- 

IS 

IS 

- 

85 

Table  4.10  FIR  Filter  Design  Parameters  and  Resulting  Filter  Order 


filter  is  extremely  high  compared  to  that  of  the  IIR  filter.  This  high  filter  order  results  in 
an  unacceptable  transient  response  as  the  filter  requires  as  many  as  85  samples  to  reach  it’s 
optimal  performance.  For  this  reason  the  FIR  filter  was  dismissed  as  a  potential  prefilter 
or  postfilter  candidate. 


4.3.4  Post-Filter  Design.  Considering  the  results  of  Section  4.3.3,  the  Post- 
Filter  was  designed  as  a  Chebyshev-I  low  pass  filter.  Since  the  Chebyshev-I  showed  only 
a  marginally  better  improvement  over  other  IIR  filter  types  for  the  prefilter,  an  extensive 
analysis  was  not  completed  for  the  postfilter.  The  following  parameters  for  the  Chebyshev- 
I  filter  were  found  to  give  the  best  results  for  the  overall  system:  The  results  for  the  overall 
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Type 

fp 

S,  (dB) 

Ss  (dB) 

order 

Data  rate 

Chebyshev-I 

80 

10 

1  Hz 

Chebyshev-I 

m 

IQ 

liQ 

80 

10 

10  Hz 

Table  4.11  HR  Post-Filter  Design  Parameters  and  Resulting  Filter  Order 
system  follow  in  Section  4.3.5 

4.3.5  Overall  Results.  In  this  section,  results  are  presented  for  the  overall  nu¬ 
merical  methods  and  digital  filtering  approach.  Figure  4.16  depicts  the  overall  approach 
involving  a  prefiltering  stage,  a  numerical  differentiator  and  a  postfiltering  stage.  The 


Figure  4.16  Overall  Block  Diagram  for  Digital  Filtering  Approach 


prefilter  and  numerical  differentiator  performance  has  been  validated  in  Sections  4.3.3  and 
4.3.2.  Based  upon  these  results,  the  overall  system  is  implemented  using  the  Chebyshev-I 
low  pass  filter  and  the  cubic  smoothing  spline  differentiator.  The  post  filter  design  param¬ 
eters  are  presented  in  Table  4.11. 


Direction 

Rate 

Mean  Error 
Colored  (m/s) 

la  (m/s) 

Mean  Error 
White  (m/s) 

la  (m/s) 

Order 

3-D  RMS 

■iia 

6.972119e-02 

6.267176e-03 

5.336845e-02 

1.080036e-02 

3-D  RMS 

6.982697e-02 

6.256767e-03 

5.354084e-02 

1.079336e-02 

3-D  RMS 

Hi 

7.859507e-02 

1.151028e-02 

7.456002e-02 

2.329147e-02 

3-D  RMS 

m 

6.841166e-02 

6.844996e-03 

6.943656e-02 

6.512700e-03 

Cubic 

3-D  RMS 

10  Hz 

6.608933e-02 

4.165794e-03 

7.414128e-02 

5.083268e-03 

1st 

3-D  RMS 

10  Hz 

6.015083e-02 

5.208232e-03 

7.061738e-02 

5.517326e-03 

2nd 

3-D  RMS 

10  Hz 

6.239051e-02 

5.323950e-03 

7.479032e-02 

6.179847e-03 

3rd 

3-D  RMS 

10  Hz 

5.262775e-02 

4.734819e-03 

6.263966e-02 

5.993183e-03 

Cubic 

Table  4.12  Overall  Performance  in  the  Presence  of  Noise  Using  Benign  Flight  Profile 
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The  overall  performance  of  the  approach  of  digital  filtering  with  numeric  differenti¬ 
ation  is  summarized  in  Table  4.12.  The  results  show  a  slight  advantage  with  the  higher 
data  rate  of  10  Hz.  The  cubic  spline  algorithm  demonstrates  the  best  overall  performance, 
including  meeting  the  specifications  of  0.005  m/s  3-D  RMS  la  for  the  case  of  colored  noise 
at  10  Hz. 

4-4  Kalman  Filtering /Smoothing 

This  section  presents  the  results  of  the  Kalman  Filtering/Smoothing  approach.  Three 
different  filter  models  are  examined:  a  constant  velocity  model,  a  constant  acceleration 
model  and  a  constant  acceleration  with  noise  modeled  as  a  first  order  Markov  process. 
Each  model  is  run  at  1  Hz  and  10  Hz.  The  benefits  of  smoothing  as  opposed  to  just 
filtering  are  presented. 

The  flight  profile  used  is  the  same  benignly  dynamic  trajectory  used  for  the  digital 
filtering  and  numeric  differentiation  approach.  The  trajectory  shows  the  initial  climb  to 
altitude  associated  with  a  takeoff.  The  ECEF  X,Y,  and  Z  velocities  are  shown  in  Figure 
4.17.  After  50  seconds,  the  aircraft  levels  off  and  assumes  a  true  constant  velocity  in  all 


Figure  4.17  X,  Y,  and,  Z  Velocities  for  Benign  Flight  Trajectory 
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directions.  Performance  is  evaluated  for  the  entire  100  second  segment  as  well  as  for  the 
last  45  seconds  which  guarantees  a  true  constant  velocity. 

4-4-1  Constant  Velocity  Model.  The  results  of  the  Constant  Velocity  Kalman 
filter /smoother  model  described  in  Section  3.2.1  are  shown  in  Tables  4.13  and  4.14.  Table 
4.13  shows  the  3-D  RMS  error  and  la  of  the  Kalman  filter /smoother  in  the  truncated 
flight  segment.  It  is  interesting  to  see  that  the  Kalman  filter  and  smoother  showed  better 
results  at  1  Hz  rather  than  10  Hz. 


Direction 

Rate 

Error  (m/s) 

la  (m/s) 

Filter  Type 

Noise  Type 

3-D  RMS 

liB 

1.513682e-02 

1.519563e-02 

Filter 

3-D  RMS 

8.156500e-03 

6.732130e-03 

Smoother 

3-D  RMS 

in 

9.072810e-03 

7.198723e-03 

Filter 

Colored 

3-D  RMS 

m 

6.706086e-03 

4.728951e-03 

Smoother 

Colored 

3-D  RMS 

Bia 

4.502794e-01 

1.226147e-01 

Filter 

White 

3-D  RMS 

2.784655e-02 

2.200576e-02 

Smoother 

White 

3-D  RMS 

10  Hz 

2.307511e-01 

3.527986e-01 

Filter 

Colored 

3-D  RMS 

10  Hz 

3.682011e-02 

2.988482e-02 

Smoother 

Colored 

Table  4.13  Constant  Velocity  Model  Using  Truncated  Flight  Profile 


When  considering  the  entire  benign  trajectory  we  must  “open  up  the  bandwidth”  of 
the  filter.  Essentially  the  non-constant  velocities  of  the  input  data  result  in  the  kinematic 
model  being  invalid.  To  tell  the  filter  to  rely  more  on  the  measurements  than  its  internal 
model,  we  increase  the  magnitude  of  the  dynamics  driving  noise,  Q.  In  this  case,  the  results 
for  the  1  Hz  and  10  Hz  data  are  very  similar.  Tuning  plots  for  the  position  and  velocity 
states  are  shown  in  Appendix  B. 

4-4-^  Constant  Acceleration  Model.  The  results  of  the  constant  acceleration 
Kalman  filter/smoother  model  described  in  Section  3.2.2  are  shown  in  Tables  4.15  and 
4.16.  Table  4.15  shows  the  3-D  RMS  error  and  Icr  of  the  Kalman  filter/smoother  in  the 
truncated  flight  segment. 

As  with  the  constant  velocity  model,  the  Kalman  filter  dynamics  driving  noise  is 
increased  to  compensate  for  the  dynamics  model  inadequacies.  These  results  are  shown  in 
Table  4.16.  Tuning  plots  for  the  position  and  velocity  states  are  shown  in  Appendix  C. 
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Direction  I  Rate 


3-D  RMS  1  Hz 
3-D  RMS  1  Hz 
3-D  RMS  1  Hz 
3-D  RMS  1  Hz 


3-D  RMS 
3-D  RMS 


Error  (m/s) 

Icr  (m/s) 

2.902978e-01 

6.008015e-02 

1.212918e-01 

2.813343e-02 

5.130680e-01 

4.753936e-02 

1.302018e-01 

2.249046e-02 

9.530139e-01 

2.477768e-02 

8.852748e-01 

3.682874e-02 

1.982058e-h00 

2.267867e-02 

1.393025e-|-00 

2.981835e-02 

Filter  Type  |  Noise  Type 


Filter 

Smoother 

Filter 

Smoother 


Filter 

Smoother 

Filter 

Smoother 


Table  4.14  Constant  Velocity  Model  Using  Entire  Benign  Flight  Profile 


Direction 

Rate 

Error  (m/s) 

1<T  (m/s) 

Filter  Type 

Noise  Type 

3-D  RMS 

1  Hz 

3.372428e-02 

2.795790e-02 

Filter 

White 

3-D  RMS 

1  Hz 

1.079198e-02 

9.758075e-03 

Smoother 

White 

3-D  RMS 

1  Hz 

2.198128e-02 

1.742112e-02 

Filter 

Colored 

3-D  RMS 

1  Hz 

1.077930e-02 

8.012129e-03 

Smoother 

Colored 

3-D  RMS 

1.041561e-02 

1.172789e-02 

Filter 

3-D  RMS 

iH 

3.921000e-03 

4.031051e-03 

Smoother 

3-D  RMS 

10  Hz 

2.585014e-02 

1.957010e-02 

Filter 

Colored 

3-D  RMS 

10  Hz 

1.389758e-02 

9.936714e-03 

Smoother 

Colored 

Table  4.15  Constant  Acceleration  Model  Using  Truncated  Flight  Profile 


Direction 

Rate 

Error  (m/s) 

la  (m/s) 

Filter  Type 

Noise  Type 

3-D  RMS 

1.618147e-01 

1.582101e-01 

Filter 

White 

3-D  RMS 

in 

4.853251e-02 

3.832201e-02 

Smoother 

White 

3-D  RMS 

HI 

7.861218e-02 

7.368700e-02 

Filter 

Colored 

3-D  RMS 

1  Hz 

2.681462e-02 

2.147360e-02 

Smoother 

Colored 

3-D  RMS 

10  Hz 

2.361016e-01 

5.188239e-01 

Filter 

White 

3-D  RMS 

10  Hz 

1.518371e-02 

1.315073e-02 

Smoother 

White 

3-D  RMS 

10  Hz 

1.180761e-01 

1.284944e-01 

Filter 

Colored 

3-D  RMS 

10  Hz 

3.151422e-02 

2.514358e-02 

Smoother 

Colored 

Table  4.16  Constant  Acceleration  Model  Using  Entire  Flight  Profile 
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4.4-3  Constant  Acceleration  with  Markov  Process  Model.  The  results  of  the 
constant  acceleration  Kalman  filter /smoother  model  with  noise  modeled  as  a  1st  order 
Markov  process  described  in  Section  3.2.3  are  shown  in  Tables  4.17  and  4.18.  Table  4.17 
shows  the  3-D  RMS  error  and  Icr  of  the  Kalman  filter/smoother  in  the  truncated  flight 
segment. 


Error  (m/s) 

la  (m/s) 

Filter  Type 

Noise  Type 

3-D  RMS 

1  Hz 

1.438570e-02 

1.443129e-02 

Filter 

White 

3-D  RMS 

1  Hz 

7.010042e-03 

6.366674e-03 

Smoother 

White 

3-D  RMS 

1  Hz 

1.051723e-02 

8.653171e-03 

Filter 

Colored 

3-D  RMS 

1  Hz 

7.704391e-03 

5.347301e-03 

Smoother 

Colored 

3-D  RMS 

10  Hz 

5.640125e-03 

8.993274e-03 

Filter 

White 

3-D  RMS 

10  Hz 

2.622584e-03 

3.179606e-03 

Smoother 

White 

3-D  RMS 

10  Hz 

1.125054e-02 

8.307786e-03 

Filter 

Colored 

3-D  RMS 

10  Hz 

9.955190e-03 

6.872084e-03 

Smoother 

Colored 

Table  4.17  Constant  Acceleration  Model  with  Markov  Process  Using  Truncated  Flight 
Profile 

As  with  the  previous  models,  the  dynamics  driving  noise,  Q,  is  increased  in  order 
to  achieve  satisfactory  performance.  Tuning  plots  for  the  position  and  velocity  states  are 
shown  in  Appendix  D. 


Direction 

Rate 

Error  (m/s) 

la  (m/s) 

Filter  Type 

Noise  Type 

3-D  RMS 

1  Hz 

1.964303e-01 

2.241990e-01 

Filter 

White 

3-D  RMS 

1  Hz 

4.448769e-02 

3.508846e-02 

Smoother 

White 

3-D  RMS 

1  Hz 

1.205297e-01 

1.247135e-01 

Filter 

Colored 

3-D  RMS 

1  Hz 

2.643877e-02 

2.122642e-02 

Smoother 

Colored 

3-D  RMS 

10  Hz 

3.652743e-01 

8.777798e-01 

Filter 

White 

3-D  RMS 

10  Hz 

1.626416e-02 

1.427469e-02 

Smoother 

White 

3-D  RMS 

10  Hz 

1.244843e-01 

1.315372e-01 

Filter 

Colored 

3-D  RMS 

10  Hz 

3.386609e-02 

2.725365e-02 

Smoother 

Colored 

Table  4.18  Constant  Acceleration  Model  with  Markov  Process  Using  Entire  Flight  Profile 
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4-5  Summary 

This  chapter  presented  the  results  of  two  distinct  approaches  to  calculating  a  refer¬ 
ence  velocity  from  discrete,  noisy  position  data:  a  digital  filtering  and  numeric  differentia¬ 
tion  approach  and  a  Kalman  filtering/smoothing  approach.  Performance  for  each  approach 
is  measured  at  data  rates  of  1  Hz  and  10  Hz,  with  a  white  noise  model  and  a  “colored” 
noise  model. 

The  results  of  these  simulations  are  encouraging.  The  performance  specification  of 
0.005  m/s  3-D  RMS  1<t  was  achieved  by  both  approaches,  although  not  for  every  case.  The 
digital  filtering  and  numeric  differentiation  method  tolerated  changes  in  dynamics  better 
than  the  Kalman  filter/smoother,  while  the  later  showed  excellent  performance  when  its 
internal  dynamics  model  was  well  matched  with  the  incoming  data. 
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V.  Conclusions  and  Recommendations 


5. 1  Overview 

This  chapter  presents  the  conclusions  derived  from  this  thesis  research  concerning  two 
different  approaches  to  determining  velocity  for  an  inverted  pseudolite-based  navigation 
reference  system:  a  digital  filtering  with  numeric  differentiation  approach  and  a  kinematic 
Kalman  filter/smoother  approach.  Also  presented  are  recommendations  for  future  research 
topics  stemming  from  this  research. 

5.2  Conclusions 

5.2.1  Digital  Filtering  and  Numerical  Methods.  In  this  research,  two  methods 
of  performing  numerical  differentiation  are  developed.  The  first  method  is  based  upon 
a  Taylor’s  series  approximation  to  the  derivative.  In  the  absence  of  noise,  the  higher 
order  differentiation  algorithms  provide  better  performance.  When  noise  is  added  to  the 
position  data,  however,  there  is  no  benefit  to  using  higher  order  numerical  differentiation 
algorithms. 

The  second  method  uses  cubic  splines  to  describe  a  curve  that  fits  the  position  data. 
The  expression  for  the  curve  is  differentiated  analytically  and  evaluated  at  the  points  of 
interest.  In  the  absence  of  noise,  this  method  performs  slightly  better  than  a  second  order 
numerical  differentiator.  In  the  presence  of  noise,  the  performance  of  the  cubic  spline  is 
comparable  to  the  Taylor’s  series-based  numerical  differentiators. 

These  numerical  differentiation  algorithms  were  tested  with  data  at  1  Hz  and  10 
Hz  rates.  Both  the  Taylor’s  series  and  cubic  spline  differentiation  algorithms  are  more 
accurate  at  the  higher  rate  of  10  Hz. 

The  data  processed  by  the  numerical  differentiation  algorithm  is  low  pass  filtered 
before  and  after  the  operation.  When  performing  a  numerical  differentiation,  it  is  beneficial 
to  remove  as  much  noise  as  possible  from  the  signal  before  differentiating.  The  prefilter, 
however,  was  not  found  to  be  particularly  effective.  Since  the  differentiation  procedure 
loses  accuracy  in  the  presence  of  noise,  improving  the  prefilter  performance  is  essential  to 
improving  the  effectiveness  of  this  approach. 
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Both  the  prefilter  and  postfilter  were  designed  as  HR  filters.  FIR  filter  designs  were 
not  able  to  meet  the  performance  requirements,  especially  in  the  area  of  allowable  passband 
loss.  The  HR  filters  were  run  forward  and  backwards  through  the  data,  allowing  the  filter 
to  achieve  a  linear  phase  response,  crucial  to  the  accurate  reconstruction  of  the  position 
data.  Given  that  the  FIR  filters  were  unable  to  meet  performance  standards  and  that  the 
HR  filters  were  non-causal,  it  is  unlikely  that  this  digital  filtering  and  numerical  methods 
approach  would  be  able  to  support  real  time  operation^. 

The  biggest  disadvantage  to  this  method  is  that  it’s  development  used  and  required 
knowledge  of  the  truth  data.  The  low  pass  filters  were  designed  to  eliminate  as  much  noise 
as  possible  without  corrupting  the  data,  and  this  was  accomplished  with  knowledge  of  the 
truth  data.  Similarly,  the  smoothing  parameter  used  in  the  cubic  smoothing  spline  algo¬ 
rithm  was  chosen  such  that  it  minimized  the  difference  between  the  estimated  velocity  and 
true  velocity.  In  a  real  world  application,  it  would  not  be  possible  to  “tune”  these  param¬ 
eters  so  exactly.  Certainly  simulation  would  be  an  aid  to  determining  these  parameters, 
however  it  remains  to  be  seen  whether  the  same  performance  levels  can  be  reached. 

In  spite  of  these  disadvantages,  this  method  was  able  to  meet  the  required  specifi¬ 
cation  of  0.005  m/s  3-D  RMS  Icr  in  simulation.  While  this  research  may  have  created 
more  questions  than  it  answered,  these  results  are  encouraging,  and  show  that  the  SARS 
concept  may  be  feasible  navigation  reference  system. 

5.2.2  Kalman  Filtering/Smoothing.  In  contrast  with  the  digital  filtering  with 
numeric  differentiation  method,  the  Kalman  filter/ smoother  did  not  demonstrate  better 
performance  at  the  higher  data  rate  of  10  Hz.  The  overall  performance  of  the  Kalman 
filter /smoother  fell  below  the  desired  accuracy  of  0.005m/s  3-D  RMS  Icr  for  a  segment  of 
the  trajectory  containing  benign  dynamics.  The  Kalman  filter /smoother  performance  was 
hindered  greatly  by  the  inadequacy  of  the  model  assumptions.  During  the  portions  of  the 
flight  profile  that  matched  the  model  assumptions,  the  Kalman  filter /smoother  performed 
well. 

^Real  time  support  is  not  required  for  the  SARS  concept.  In  this  post  processing  environment  it  is 
prudent  to  take  advantage  of  non-causal  filters. 
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One  way  to  improve  the  performance  of  Kalman  filter  would  be  to  lower  the  mea¬ 
surement  noise.  Since  the  position  accuracy  of  0.1  m  is  given  in  the  problem  statement, 
it  may  not  be  directly  possible  to  lower  the  measurement  noise.  However,  one  could  alter 
the  Kalman  filter  models  to  accept  ranging  measurements.  In  this  manner,  ranging  error 
terms  such  as  tropospheric  delay  and  multipath  could  be  estimated  and  potentially  re¬ 
duced.  Since  the  ranging  equations  are  non-linear,  an  extended  Kalman  filter  formulation 
or  other  compensation  may  be  required. 

To  exploit  the  availability  of  past  and  future  measurements  in  the  Kalman  filter, 
a  fixed  interval  Kalman  smoother  was  implemented.  While  there  was  a  benefit  to  the 
smoothing  operation,  it  was  not  enough  to  completely  overcome  the  inherent  model  inad¬ 
equacy.  The  Fraser  smoothability  criteria  indicates  that  the  kinematic  models  used  in  this 
research  are  not  very  smoothable  when  the  dynamics  driving  noise  is  low  (i.e.  the  internal 
dynamics  model  matches  the  observations)  and  this  is  reflected  in  the  results.  When  Q  is 
increased  to  compensate  for  an  inadequate  model,  the  benefit  of  smoothing  is  greater. 

While  the  performance  objectives  were  not  met  using  the  Kalman  filter /smoother 
approach  for  the  entire  benign  trajectory,  there  are  advantages  to  this  method.  Foremost 
is  that  this  method  is  less  reliant  upon  knowledge  of  the  truth  data  in  order  to  “tune” 
the  system,  as  long  as  one  has  confidence  in  the  dynamics  model.  Since  the  Kalman 
filter/smoother  showed  promise  when  the  model  was  well  matched  to  the  measurement 
data,  it  still  may  be  a  useful  research  avenue  for  SARS. 

5.3  Recommendations 

The  following  are  just  a  few  of  the  many  possible  research  topics  for  future  considera¬ 
tion.  The  recommendations  presented  represent  research  areas  most  important  to  proving 
the  feasibility  of  the  SARS  concept. 

5.3.1  Receiver  Geometry.  As  was  exhibited  in  this  research,  proper  receiver 
location  is  essential  to  maintaining  good  geometry  (e.g.  DOP)  over  the  course  of  a  flight. 
Since  the  accuracy  of  the  solution  to  the  GPS  ranging  equations  worsens  as  the  receiver 
locations  become  more  coplanar,  it  is  likely  that  there  exists  a  limit  on  what  the  acceptable 
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DOP  is  to  maintain  a  position  accuracy  of  0.1  m.  If  this  DOP  limit  can  be  found,  it  may 
be  possible  to  determine  an  “optimal”  placement  of  pseudolite  receivers  that  maintains 
good  geometry  throughout  a  flight  profile. 

5.3.2  Flight  Profiles.  Flight  profiles  used  to  test  navigation  systems  are  typically 
2-3  hour  missions  involving  a  combination  of  benign  and  high  dynamic  segments  [39].  This 
research  also  showed  that  such  flight  profiles  will  have  to  be  modified  to  maintain  good 
receiver  geometry.  It  may  be  necessary  to  consider  the  problem  of  determining  appropriate 
flight  profiles  and  optimizing  receiver  placement  simultaneously. 

5.3.3  Error  Modeling.  For  the  unique  inverted  pseudolite  concept  of  SARS,  creat¬ 
ing  SARS-specific  tropospheric  and  multipath  models  could  increase  positioning  accuracy, 
and  subsequently  increase  velocity  accuracy.  Using  empirical  data  it  may  be  possible  to 
quantify  and  model  the  tropospheric  delay  and  multipath  error  terms,  given  the  relatively 
consistent  climate  of  the  New  Mexico  desert,  as  well  as  apriori  knowing  the  receiver  loca¬ 
tions. 

5.3.4  Adaptive  Kalman  Filtering.  This  research  shows  the  potential  for  a  kine¬ 
matic  Kalman  filter  based  solution  if  the  underlying  filter  model  adequately  describes  the 
dynamics  of  the  trajectory.  While  one  Kalman  filter  may  not  be  able  to  meet  the  accuracy 
needed  by  SARS,  it  may  be  worthwhile  to  consider  an  adaptive  filtering  technique.  For 
example,  a  Multiple  Model  Adaptive  Estimation  (MMAE)  [25]  technique  could  be  used  to 
set  up  a  bank  of  Kalman  filters  which  run  in  parallel.  Using  residual  monitoring  techniques 
the  outputs  of  these  Kalman  filters  can  be  multiplexed  and/or  blended  in  an  attempt  to 
achieve  better  performance  than  that  of  a  single  Kalman  filter. 
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Appendix  A.  Calculation  of  RMS  Error 


The  objective  of  this  appendix  is  to  illustrate  the  error  metrics  used  in  this  thesis 
research.  A  scalar  metric,  the  Root  Mean  Square  (RMS),  is  used  to  describe  the  error  of  a 
vector  of  computed  values  from  the  true  values.  Typically  these  vectors  are  a  time  history 
of  position  or  velocity. 

Consider  a  time  history  of  three  dimensional  position  expressed  as  three  vectors,  X, 
Y  and  Z,  as  well  as  the  estimates  of  position  X,  Y  and  Z,  all  of  length  k.  The  3-D  RMS 
error  metric  is  given  by: 

i  ^  ^J{Xiu)  -  Xit,)y  +  (y(t,)  -  Y{t,)y  +  (z{t,)  -  z{u)y  (a.i) 

i  =  l 


In  the  case  of  data  from  more  than  one  simulation  (Monte  Carlo  analysis),  it  is  also 
useful  to  compute  the  following  mean  and  standard  deviation  statistics: 


mRMS 


^'^mRMS 


i.f;RMS(i) 


i=l 


\  ^•E(RMS(0-mRMS)2 

\  ”  «=i 


(A.2) 

(A.3) 


where  mRMS  denotes  the  mean  RMS  error,  lOmRMS  denotes  the  standard  deviation  in 
RMS  error  and  n  represents  the  number  of  simulation  runs  performed. 
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Appendix  B.  Constant  Velocity  Kalman  Filter  Results 

This  appendix  contains  the  tuning  plots  for  the  Constant  Velocity  Kalman  filter.  A 
legend  for  the  plots  is  presented  below: 


Table  B.l  Legend  for  Filter  Tuning  Plots 


Symbol 

Definition 

-  Solid  Line 

Mean  Error 

•  •  •  Dotted  Line 

Mean  Error  ±  True  Sigma 

—  Dashed  Line 

±  Filter  Predicted  Sigma 

X  Position  Tuning  Y  Position  Tuning  Z  Position  Tuning 


ime  (sec)  Time  (sec)  Time  (sec) 


X  Position  Tuning  Y  Position  Tuning  Z  Position  Tuning 


Figure  B.2  White  Noise  at  1  Hz,  Proper  Model 
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ime  (sec)  Time  (sec)  Time  (sec) 


X  Position  Tuning  Y  Position  Tuning  Z  Position  Tuning 


Figure  B.3  Colored  Noise  at  10  Hz,  Proper  Model 
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Figure  B.4  White  Noise  at  10  Hz,  Proper  Model 
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Figure  B.6  White  Noise  at  1  Hz,  Improper  Model 
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Appendix  C.  Constant  Acceleration  Kalman  Filter  Results 

This  appendix  contains  the  tuning  plots  for  the  Constant  Velocity  Kalman  filter.  A 
legend  for  the  plots  is  presented  below: 


Table  C.l  Legend  for  Filter  Tuning  Plots 


X  Position  T uning  Y  Position  T uning  Z  Position  Tuning 
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Kalman  Filter 
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Appendix  D.  Constant  Acceleration,  1st  Order  Markov  Kalman  Filter  Results 


This  appendix  contains  the  tuning  plots  for  the  Constant  Velocity  KaJman  filter.  .4 
legend  for  the  plots  is  presented  below; 

Table  D.l  Legend  for  Filter  Tuning  Plots 


Symbol 

Definition 

-  Solid  Line 

Mean  Error 

Mean  Error  x  True  Sigma 

—  Dashed  Line 

±  Filter  Predicted  Sigma 
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